A reservoir operation state recognition method and device coupling a hidden Markov model and a decision tree

Abstract

The application discloses a reservoir operation state recognition method and device coupling a hidden Markov model and a decision tree, and the method comprises the following steps: obtaining reservoir operation observation data of a target reservoir, and constructing a time sequence containing multi-dimensional hydrological variables; a hidden Markov model is constructed, parameter estimation of the hidden Markov model is carried out based on the time sequence, the model after the parameter estimation is used to decode the time sequence, and a hidden state sequence representing a reservoir operation mode is obtained; the hidden state sequence is taken as a label, the time sequence is taken as a feature, a decision tree model is constructed and trained, and reservoir operation discrimination rules corresponding to different hidden states are extracted; the optimal hidden state number is determined based on preset model optimization indexes, and the hidden state sequence corresponding to the optimal hidden state number and the reservoir operation discrimination rules are output. The application solves the problem that an existing model cannot accurately capture the periodicity of reservoir scheduling, and realizes fine identification of the operation state and analysis of physical rules.

Description

Technical Field

[0001] This invention belongs to the field of water conservancy engineering, and in particular, it relates to a method and device for identifying the operation status of a reservoir by coupling a hidden Markov model and a decision tree. Background Technology

[0002] Reservoir operation is a complex and dynamic process influenced by hydrological and meteorological conditions, engineering physical constraints, and human-made scheduling rules. Accurately identifying the operational status of reservoirs from long-term hydrological observation data is the technical foundation for inverting scheduling rules, assessing scheduling compliance, and analyzing downstream hydrological responses. By mining potential patterns in time-series data such as water level and flow, continuous hydrological processes can be discretized into scheduling stages with specific physical meanings, providing quantitative basis for refined watershed management.

[0003] Existing methods for reservoir operation status identification mainly include rule-matching methods based on scheduling charts, statistical methods based on cluster analysis, and data-driven methods based on machine learning. Among these, Hidden Markov Models (HMMs) are widely used in reservoir status identification because they can infer unobservable potential operation modes from observable hydrological variables. These methods typically assume that the reservoir operation status follows a Markov chain, describe the switching probabilities between different scheduling stages (such as water storage and flood discharge) through a state transition matrix, and use the probability distribution of observed variables to characterize the hydrological features of each state.

[0004] However, existing technologies have limitations in handling reservoir scheduling problems with strong seasonal characteristics, mainly manifested in the mismatch between model assumptions and physical laws, and the lack of physical meaning in state determination. Specifically, standard Hidden Markov Models (HMMs) typically assume that state transition probabilities are time-invariant (homogeneity assumption), meaning that the probability of a reservoir switching states at any time of year is constant. This ignores the strong periodicity of reservoir scheduling driven by flood season water level constraints and inflow conditions, making it difficult for the model to capture the scheduling logic that dynamically changes with seasons and operating conditions. In addition, existing methods often rely on purely statistical criteria such as the Akaike Information Criterion (AIC) or the Bayesian Information Criterion (BIC) to determine the number of hidden states, without considering the dynamic differences between states. This can easily lead to the model identifying statistically separable but physically similar redundant states, reducing the interpretability and physical rationality of the identification results. Summary of the Invention

[0005] The purpose of this invention is to provide a method and apparatus for identifying the operational status of a reservoir by coupling a hidden Markov model and a decision tree, so as to solve the above-mentioned problems existing in the prior art.

[0006] Technical solution: A method for identifying the operational status of a reservoir by coupling a hidden Markov model and a decision tree, comprising:

[0007] Obtain reservoir operation observation data of the target reservoir and construct a time series containing multidimensional hydrological variables;

[0008] A hidden Markov model is constructed, and the parameters of the hidden Markov model are estimated based on the time series. The time series is then decoded using the model with the estimated parameters to obtain the hidden state sequence that represents the reservoir operation mode.

[0009] Using the hidden state sequence as the label and the time series as the feature, a decision tree model is constructed and trained to extract the reservoir operation discrimination rules corresponding to different hidden states.

[0010] The optimal number of hidden states is determined based on the preset model optimization index, and the hidden state sequence corresponding to the optimal number of hidden states and the reservoir operation discrimination rules are output.

[0011] According to one aspect of this application, a reservoir operation status identification device coupled with a hidden Markov model and a decision tree includes:

[0012] Memory, used to store computer programs;

[0013] The processor is used to execute a computer program stored in the memory, and when the computer program is executed, it implements the reservoir operation status identification method that couples the hidden Markov model and the decision tree as described in the above technical solution.

[0014] Beneficial effects: This invention replaces the traditional constant matrix with a time-varying transition probability matrix driven by covariates, solving the problem that existing models cannot accurately capture the periodic patterns of reservoir scheduling, and realizing refined identification of operating status and analysis of physical rules. Attached Figure Description

[0015] Figure 1 A flowchart illustrating the steps of a reservoir operation status identification method that couples a hidden Markov model and a decision tree, as provided in an embodiment of this application.

[0016] Figure 2 This is a schematic diagram of the Weibull distribution fitting results of the observed variables in the Hidden Markov Model provided in the embodiments of this application.

[0017] Figure 3 This is a schematic diagram of the Gaussian distribution fitting results of the observed variables in the Hidden Markov Model provided in the embodiments of this application.

[0018] Figure 4 A flowchart illustrating the steps for extracting reservoir operation discrimination rules corresponding to different hidden states, as provided in this application embodiment.

[0019] Figure 5 A flowchart illustrating the steps for performing input variable screening based on variable importance feedback, as provided in an embodiment of this application.

[0020] Figure 6 This is a schematic diagram showing the ranking results of different input combination variables provided in the embodiments of this application.

[0021] Figure 7 A flowchart illustrating the steps of solving the regression coefficients of a multinomial Logistic regression model using the extended expectation maximization algorithm, as provided in this application embodiment. Detailed Implementation

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

[0023] It should be noted that the terms include and have, and any variations thereof, are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or device that includes a series of steps or units is not necessarily limited to those steps or units that are explicitly listed, but may include other steps or units that are not explicitly listed or that are inherent to such process, method, product, or device.

[0024] like Figure 1 As shown, a method for identifying the operational status of a reservoir by coupling a hidden Markov model and a decision tree includes the following steps:

[0025] Obtain reservoir operation observation data for the target reservoir and construct a time series containing multidimensional hydrological variables.

[0026] In this embodiment, the reservoir operation observation data mainly consists of a set of physical quantities reflecting the daily operation status of the reservoir and the hydrological situation of the watershed. Specifically, the acquired data typically includes, but is not limited to, the inflow rate Q. in Outbound flow Q out The data includes the reservoir water level (WL) and the reservoir storage volume (Volume). The process of constructing a time series generally involves organizing, cleaning, and aligning the scattered observation data in chronological order to form a multidimensional matrix. Each row represents a time step, and each column represents a hydrological variable. Constructing a multidimensional hydrological variable time series is the foundation for subsequent model training. In some embodiments, flow data from downstream control stations can be introduced as auxiliary variables to support data cleaning. For example, for Reservoir A in the Yangtze River basin, daily hydrological data from 2015 to 2023 can be collected to form a long-sequence sample set.

[0027] A hidden Markov model is constructed, and parameters of the hidden Markov model are estimated based on time series data. The time series data is then decoded using the model with the estimated parameters to obtain the hidden state sequence representing the reservoir operation mode.

[0028] Specifically, Hidden Markov Models (HMMs) are used to describe Markov processes containing hidden unknown parameters. In the context of reservoir operation, the explicit observed variables are flow rate and water level, while the reservoir's operational mode is a hidden state that cannot be directly observed. These operational modes include flood control, water supply, and ecological operation modes. Constructing an HMM establishes a probabilistic relationship between observed variables and hidden states. Parameter estimation involves using the Expectation-Maximization (EM) algorithm or its variants to infer parameters such as state transition probabilities and emission probabilities from observed data, maximizing the likelihood probability of the model generating the observed sequence. Decoding involves using the Viterbi algorithm or other optimal path search algorithms to infer the hidden state sequence most likely to produce the observed sequence, given the model parameters and the observed sequence. The resulting hidden state sequence discretizes the continuous time axis into several state intervals with specific physical meanings, revealing the temporal evolution of reservoir operation.

[0029] Using the hidden state sequence as labels and the time series as features, a decision tree model is constructed and trained to extract the reservoir operation discrimination rules corresponding to different hidden states.

[0030] In this embodiment, since the hidden states identified by HMM are usually represented by numerical labels, which lack intuitive physical interpretability, a decision tree model, such as Classification and Regression Tree (CART), can be introduced for rule extraction.

[0031] Specifically, the hidden state sequence is used as the target variable, and the corresponding hydrological variable (Q) is used as the target variable. in Q out A supervised learning dataset is constructed using features (WL, Volume) as feature variables. By training a decision tree, the key features and threshold conditions that determine the occurrence of a certain hidden state can be extracted. The extracted discrimination rules are usually expressed in the form of if-then logic. For example, if the water level is higher than 162 meters and the inflow is less than 1000 cubic meters per second, then the reservoir is in a water supply scheduling state. The rule-based output allows schedulers to intuitively understand the hydrological boundary conditions and scheduling logic represented by each hidden state.

[0032] The optimal number of hidden states is determined based on the preset model optimization index, and the hidden state sequence corresponding to the optimal number of hidden states and the reservoir operation discrimination rules are output.

[0033] Specifically, the number of hidden states in a Hidden Markov Model (HMM), i.e., the K value, is a hyperparameter that needs to be optimized based on data characteristics and application requirements. Pre-defined model optimization metrics can be statistical indicators, such as AIC and BIC, or custom metrics incorporating physical meaning, such as state redundancy. Determining the optimal number of hidden states typically involves modeling and evaluating different K values, such as K=2, 3, 4, and 5, and selecting the model with the optimal metrics as the final result. The output includes not only the hidden state sequence under the optimal K value but also the decision tree discrimination rules trained based on that K value. This ensures that the final identified reservoir operation state is both statistically significant and conforms to the physical laws of actual scheduling, avoiding overly coarse or fine state divisions.

[0034] In other words, a reservoir operation status identification method that couples a hidden Markov model with a decision tree can also be:

[0035] Acquire reservoir operation observation data for the target reservoir and construct a time series containing multidimensional hydrological variables. For a predetermined number of candidate hidden states, construct a Hidden Markov Model (HMM) for each. Estimate the parameters of the HMM based on the time series and decode the time series using the estimated model to obtain the hidden state sequence representing the reservoir operation mode corresponding to each candidate hidden state number. Use the hidden state sequence corresponding to each candidate hidden state number as a label and the time series as a feature to construct and train decision tree models to extract reservoir operation discrimination rules corresponding to different hidden states. Determine the optimal hidden state number from the predetermined number of candidate hidden states based on a preset model optimization index, and output the hidden state sequence corresponding to the optimal hidden state number and the reservoir operation discrimination rule.

[0036] One possible implementation involves constructing a time series containing multidimensional hydrological variables, and also includes imputing missing values ​​and filtering the series data of reservoir operation observations, specifically:

[0037] For outflow data with missing values ​​in reservoir operation observation data, the same-period flow data of the nearest downstream river hydrological station is obtained, and a linear regression model of reservoir outflow and river hydrological station flow for the period without missing values ​​is established. The missing values ​​are then calculated and filled based on the linear regression model.

[0038] In this embodiment, reservoir outflow data may be missing due to equipment failure or transmission interruption. Directly removing periods containing missing values ​​may result in information loss, especially in continuous modeling. Therefore, a correlation-based interpolation imputation method is adopted.

[0039] Specifically, for Reservoir A, the nearest downstream hydrological station B is selected as the reference station. The outflow Q from Reservoir A during historical periods with no missing data is used. out With the flow rate Q at hydrological station B stationPerform regression analysis to establish a linear regression model. For example, the linear regression equation obtained after fitting could be:

[0040] Q out =1.02*Q station +263;

[0041] The coefficients 1.02 and intercept 263 are regression parameters estimated using the least squares method.

[0042] Based on a linear regression model, when the outflow from Reservoir A is missing but the flow at Hydrological Station B is known, the missing value can be calculated by substituting it into the above formula. This embodiment utilizes the strong physical correlation between upstream and downstream flows, and compared to simple mean filling or linear interpolation, it has higher physical reliability.

[0043] The length of the imputed time series is validated. If the length of a continuous time series is less than the preset minimum sample length threshold, the sequence is discarded, and only the time series whose length meets the modeling requirements are retained for constructing the Hidden Markov Model.

[0044] In this embodiment, parameter estimation of the Hidden Markov Model relies on a sufficiently long time series to capture state transition patterns. Too short a sequence segment may lead to unstable or non-convergent estimation of the state transition probability matrix. Therefore, a minimum sample length threshold is set for filtering.

[0045] For example, the minimum sample length threshold can be preset to 100 days. After missing value imputation, each continuous time series is examined. If the length of a continuous data segment is less than 100 days, it is considered an invalid segment and discarded. Only continuous sequences with a length greater than or equal to 100 days are retained for subsequent HMM training. This ensures that the data input to the model has a sufficient time span to encompass the complete hydrological event process, improving the robustness of the model.

[0046] In an exemplary embodiment, constructing a hidden Markov model includes constructing a response model, specifically:

[0047] Identify the probability density distribution characteristics of multidimensional hydrological variables, and construct response models using non-normal distribution functions for variables with extreme value distributions.

[0048] Alternatively, it can be described as identifying the probability density distribution characteristics of multidimensional hydrological variables, and using non-normal distribution functions to construct response models for variables with long-tail characteristics.

[0049] Specifically, the response model of a Hidden Model (HMM), i.e., the emission probability model, determines the distribution of observed data under a given latent state. Traditional methods often assume that observed variables follow a Gaussian distribution, i.e., a normal distribution. However, hydrological data typically exhibit significant skewness and long-tail characteristics, meaning that extreme events, though small in probability, have a large impact. By performing distribution fit tests on variables such as inflow and outflow rates and water levels, for example using the AIC criterion or the Kolmogorov-Smirnov (KS) test, non-normal characteristics can be identified. To address non-normal characteristics, the assumption of a single Gaussian distribution is abandoned, and a non-normal distribution function that better conforms to hydrophysical laws is adopted to construct the response model. This requires extending the HMM modeling framework programmatically to define new distribution classes to support parameter estimation and probability calculation for non-normal distributions.

[0050] Among them, the non-normal distribution function includes the Weibull distribution function and the log-normal distribution function. The Weibull distribution function is used to describe the probability distribution of the inflow variable, and the log-normal distribution function is used to describe the probability distribution of the outflow, reservoir water level and reservoir storage variables.

[0051] In this embodiment, specific distribution types are specified based on the physical properties and statistical characteristics of hydrological variables. Specifically, the inflow rate Q... in Due to the suddenness and long tail of rainfall-runoff processes, the Weibull distribution is chosen for description. The probability density function of the Weibull distribution can flexibly adapt to various forms from exponential to normal distributions, making it particularly suitable for describing extreme flow rates. The outflow flow Q... out The reservoir water level (WL) and storage capacity (Volume) are strongly constrained by scheduling rules, and their values ​​typically fluctuate within a certain range and must be positive. Therefore, a log-normal distribution is chosen for description. The log-normal distribution, by taking the logarithm of the variables, follows a normal distribution, effectively handling data skewness. Through targeted distributional assumptions, the model's likelihood estimation of the observed data becomes more accurate, improving the precision of state identification. For example... Figure 2 and Figure 3 As shown, taking the data from this embodiment as an example, the inbound flow Q is compared between the Weibull distribution (a newly defined distribution type) and the Gaussian distribution embedded in the model. in The results showed that the former's fitting results were significantly better than the latter's.

[0052] Furthermore, when constructing the response model, the shape and scale parameters of the Weibull distribution function and the log-normal distribution function are initialized using the maximum likelihood estimation method and the moment estimation method, respectively.

[0053] In this embodiment, the parameter estimation of the Hidden Markov Model (HMM) typically employs the EM algorithm, which is sensitive to initial parameters. To accelerate convergence and avoid getting trapped in local optima, the parameters of the response model need to be properly initialized. For the Weibull distribution, it includes shape and scale parameters. This embodiment uses the Maximum Likelihood Estimation (MLE) method to estimate the initial values ​​of these two parameters based on the observed data.

[0054] Specifically, a log-likelihood function for the Weibull distribution is constructed, and the parameter value that maximizes the likelihood function is used as the initial value through a numerical optimization algorithm. For the log-normal distribution, which includes the logarithmic mean and logarithmic standard deviation, the method of moments is employed. Using the sample mean and sample variance of the observed data, the initial values ​​of the distribution parameters are directly solved using the relationship between moments. This statistically based initialization strategy provides a good starting point for subsequent complex HMM iterative training. In the specific code implementation, the initialization of the above parameters and subsequent iterative optimization process can be completed by calling the optimization function in conjunction with the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm.

[0055] This embodiment can solve the common problem of missing values ​​in hydrological data, as well as the problem that the conventional normal distribution assumption cannot adapt to the long-tail characteristics of hydrological data, providing a high-quality data foundation and accurate emission probability description data for subsequent HMM modeling.

[0056] In a preferred implementation, the Hidden Markov Model is a non-homogeneous Hidden Markov Model, which contains a time-varying state transition probability matrix, and the state transition probabilities in the time-varying state transition probability matrix are functions of time covariates.

[0057] The time covariates should include at least the inflow factor Q, which reflects the characteristics of the watershed's inflow. in (t), the outflow factor Q, which reflects the intensity of reservoir discharge. out (t), water level factor WL(t) reflecting the current operating potential energy of the reservoir, and reservoir capacity factor Volume(t) reflecting the reservoir capacity occupancy status;

[0058] By establishing a dynamic mapping relationship between state transition probabilities and time covariates, the process of identifying hidden state sequences can be adapted to the seasonal and periodic characteristics of reservoir scheduling.

[0059] In other words, the Hidden Markov Model (HMM) is a non-homogeneous HMM. The HMM includes a time-varying state transition probability matrix, where the state transition probabilities are not fixed constants but functions of time covariates. When constructing the HMM, physical factors reflecting the reservoir's hydrological conditions and operational status are introduced as time covariates. These time covariates include at least the inflow factor Q, which reflects the characteristics of the watershed's inflow. in(t), the outflow factor Q, which reflects the intensity of reservoir discharge. out The water level factor WL(t) reflects the current operating potential energy of the reservoir, and the reservoir capacity factor Volume(t) reflects the occupancy status of the reservoir.

[0060] In this embodiment, the traditional homogeneous HMM assumes that the transition probability matrix A is constant, meaning that the probability of a reservoir transitioning from a flood discharge state to a water storage state is the same regardless of whether it is the flood season or the dry season, which obviously violates the laws of physics. Therefore, a non-homogeneous HMM is constructed, in which the state transition probability matrix A(t) changes dynamically with time t.

[0061] Specifically, define the time covariate vector X at time t. t for:

[0062] X t =[1, Q in (t), Q out [(t), WL(t), Volume(t)];

[0063] Where 1 is the constant corresponding to the intercept term, Q in Q(t) represents the inflow rate at time t. out WL(t) is the outflow rate at time t, WL(t) is the water level at time t, and Volume(t) is the water storage capacity at time t.

[0064] These physical factors, as driving variables, directly determine the tendency of the reservoir's operating state to shift in the next moment. For example, when the water level WL(t) is close to the flood control limit and the inflow Q... in As (t) increases, the probability of the model predicting a shift to flood discharge mode will increase, demonstrating its adaptability to changes in operating conditions.

[0065] In a further embodiment, a dynamic mapping relationship is constructed using a multinomial logistic regression model. The specific calculation process includes:

[0066] Define the time covariate vector X at time t. t For Q in (t), Q out Combinations of (t), WL(t), and Volume(t);

[0067] Calculate the linear predictor η for the transition from hidden state i to hidden state j at time t. ij (t), the calculation formula is:

[0068] ;

[0069] Where, β ij，0 For the baseline intercept term, β ij，1 To βij，4 These are the regression coefficients for the inflow factor, outflow factor, water level factor, and reservoir capacity factor, respectively; the regression coefficients are determined by parameter estimation of the hidden Markov model.

[0070] Based on linear predictor η ij (t), calculate the state transition probability P(S) from hidden state i to hidden state j at time t. t =j|S t-1 =i,X t The calculation formula is:

[0071] ;

[0072] Among them, S t S represents the hidden state at time t. t-1 Let K represent the hidden state at time t-1, and K be the total number of hidden states. x It is an exponential function, and k is the summation index variable.

[0073] In this embodiment, the continuously changing covariates are mapped to the probability space of the interval [0, 1] by multinomial logistic regression, and the sum of the transition probabilities of all target states is guaranteed to be 1.

[0074] Specifically, for each initial state i, K-1 regression models need to be constructed to predict the tendency to transition to different target states j. Linear predictor η ij (t) is essentially a linear combination of covariates, and its magnitude reflects the relative potential energy at which the transition occurs. The regression coefficient β ij These are the core parameters that the model needs to learn, quantifying the degree of influence of each physical factor on the state transition. For example, if β ij，3 A positive value for the coefficient of the water level factor indicates that the higher the water level, the greater the probability of transitioning from state i to state j.

[0075] In a further embodiment, when constructing a multinomial logistic regression model, to ensure the identifiability of the model parameters, for each initial hidden state i, a reference target hidden state k=1 is set, and the linear predictor corresponding to the reference target hidden state is constantly set to 0, i.e., η is set. i1 (t)=0; the linear predictor η of the remaining target hidden state j (j≠1) ij (t) represents the logarithmic occurrence ratio relative to the hidden state of the reference target.

[0076] Specifically, due to the limitations of probability normalization, estimating the parameters independently for all K target states would lead to an infinite number of solutions to the equations, i.e., over-parameterization, making the model unidentifiable. Therefore, a baseline category logic strategy is adopted. Specifically, for any initial state i, target state 1 is manually designated as the reference state, and its corresponding regression coefficient vector is forced to be all 0, i.e., the corresponding linear predictor η... i1 (t) is always equal to 0. At this time, the linear predictor η of other target states j... ij (t) actually represents the natural logarithm of the ratio of the probability of transitioning to state j to the probability of transitioning to state 1, i.e., the log-odds ratio. This eliminates redundant degrees of freedom and ensures the model parameter β... ij The unique solution guarantees the stability of numerical calculations.

[0077] like Figure 7 As shown, according to one aspect of this application, parameter estimation of a hidden Markov model is performed, specifically by using the extended expectation-maximization algorithm to solve for the regression coefficients of a multinomial logistic regression model, including:

[0078] In the expectation step, based on the time series, the posterior probability of the hidden state is calculated using the forward-backward algorithm, where the state transition probability in the calculation process adopts the time-varying value obtained by the linear predictor at the current time.

[0079] In the maximization step, a weighted log-likelihood function is constructed based on the posterior probability, and the regression coefficients are iteratively updated using the iterative weighted least squares method or the quasi-Newton method until the model converges.

[0080] In this embodiment, since the transition probability changes over time, the standard Baum-Welch algorithm is no longer applicable, and the extended EM algorithm must be used. In the expectation step, the forward variable α is calculated. t (j) and backward variable β t In case (i), the transition probability used in the recursive formula is no longer a constant a. ij It is not the specific time-varying value a ij (t). For example, the recursive formula for the forward variable is revised as follows:

[0081] ;

[0082] Where b j (O t ) is the emission probability.

[0083] In the maximization step, the goal is to update the regression coefficient β, which is equivalent to solving a weighted logistic regression problem, where the sample weights are the expected number of transitions ξ calculated in the expectation step. t(i, j). Since this optimization problem has no closed-form solution, iterative weighted least squares or quasi-Newton methods are preferred for numerical solution. Multiple iterations are performed until the increment of the log-likelihood function is less than a preset threshold, such as 1e. -6 .

[0084] In one embodiment of this application, the time series is decoded using the model after parameter estimation. Specifically, the Viterbi algorithm is used to obtain the optimal hidden state sequence. In the recursive process of the Viterbi algorithm, the state transition probability at time t directly calls the element value in the time-varying state transition probability matrix A(t) calculated at that time.

[0085] Specifically, the goal of decoding is to find the hidden state path with the highest probability. The standard Viterbi algorithm uses a fixed A matrix during dynamic programming optimization. In this embodiment, the Viterbi algorithm is modified to a time-varying Viterbi algorithm. The maximum path probability δ to state j at time t is calculated. t When (j), the transition probability used depends on the covariate X at the current time step. t Calculated a ij (t). The recursive formula is:

[0086] δ t (j)=max (1≤i≤K) [δ t-1 (i)·a ij (t)]·b j (O t );

[0087] Among them, b j (O t O is the observed value at time t in hidden state j. t The probability of emission.

[0088] In this way, the decoding process fully considers the hydrological boundary conditions at the time. For example, during high water levels in the flood season, a ij (t) will automatically reduce the probability of transitioning to the water storage state, making the decoded state sequence more consistent with the actual logic of flood control scheduling.

[0089] This embodiment addresses the shortcomings of traditional models in adapting to seasonal changes in reservoir scheduling by introducing physical covariates to drive the time-varying evolution of state transition probabilities.

[0090] like Figure 4 As shown, in one possible embodiment, a decision tree model is constructed and trained to extract reservoir operation discrimination rules corresponding to different hidden states, specifically including:

[0091] The state labels in the hidden state sequence are merged with the corresponding multidimensional hydrological variables and time covariates to construct a supervised learning sample set.

[0092] Based on the pre-defined reservoir scheduling cycle characteristics, the supervised learning sample set is divided into a flood season sample set and a non-flood season sample set;

[0093] Pre-configured classification and regression tree models were trained for flood season and non-flood season sample sets respectively. The Gini index was used as the splitting criterion to extract the reservoir operation discrimination rules for different time periods.

[0094] In this embodiment, reservoir scheduling exhibits seasonal variations. For example, during the flood season, typically from June to September, the primary objective is flood control, with strict water level control; while during the non-flood season, from October to May of the following year, the main objectives are water storage and supply, with water levels maintained at as high a level as possible. Training a single decision tree using all the data could lead to rule conflicts or ambiguity. Therefore, a divide-and-conquer strategy is preferred.

[0095] Specifically, construct a structure containing [State, Q] in Q out The total sample set is defined as [WL, Volume], where State is the hidden state label of the reservoir operation output by the HMM. Based on the month attribute of the date, the data is physically split into two subsets: the flood season sample set and the non-flood season sample set.

[0096] Furthermore, a CART decision tree is constructed for each subset. During training, the optimal splitting feature and splitting threshold are selected using the Gini Index minimization principle. For example, a decision tree trained on a flood season sample set might generate the following rule: if the water level WL > 158.5m and the inflow Q... in >5000m 3 If the flow rate is / s, it is classified as State 3, i.e., high-flow-rate flood discharge mode. For the non-flood season sample set, the rule might be: if the water level WL ≥ 162m and the outflow Q... out <1000m 3 If the value is / s, it is determined to be State 1, i.e., water storage and conservation mode, which can also be a low-flow operation mode. By extracting rules in different time periods, the resulting rule set is not only more accurate, but also conforms to the phased scheduling concept in reservoir scheduling procedures, making it easier for scheduling personnel to understand and verify. In some optional implementation methods, more time periods, such as water supply period and ecological scheduling period, can be further subdivided according to the specific reservoir scheduling procedures for separate modeling.

[0097] This embodiment solves the problem that a single rule model cannot simultaneously cover the complex scheduling logic of flood season and non-flood season by dividing physical time periods, thus improving the interpretability of the model.

[0098] In one embodiment of this application, the optimal number of hidden states is determined based on a preset model optimization index, including calculating the inter-state transition probability difference, which is used to measure the distinguishability of any two hidden states in their transition behavior.

[0099] Calculate the difference in state transition probabilities D between any two hidden states i and j. trans (i, j), the calculation formula is:

[0100] ;

[0101] Where K is the total number of hidden states, k is the summation index variable, and a ik Let a represent the probability of transitioning from hidden state i to hidden state k. jk Let a represent the probability of transitioning from hidden state j to hidden state k. ki Let a represent the probability of transitioning from hidden state k to hidden state i. kj Let represent the probability of transitioning from hidden state k to hidden state j, and || denotes the absolute value operation. The first term of the formula measures the difference in transition behavior between hidden states i and j, and the second term measures the difference in transition behavior between hidden states i and j.

[0102] In this embodiment, the difference in state transition probabilities D trans It considers not only the differences in the transition behaviors between the two states, but also the differences in the transition behaviors between the two states.

[0103] Specifically, the first term in the formula for calculating the difference in transition probabilities between states actually calculates half the Manhattan distance between the vectors in the i-th and j-th rows of the transition matrix, reflecting the similarity in the transition tendencies of the two states at the next time step. The second term in the formula calculates the difference between the vectors in the i-th and j-th columns of the transition matrix, reflecting the difference in the ease with which the system enters these two states.

[0104] In another embodiment of this application, the optimal number of hidden states is determined based on a preset model optimization index, including calculating the difference in transition probabilities between states:

[0105] The time-varying state transition probability matrix is ​​averaged over time to obtain the average state transition probability matrix. Based on the average state transition probability matrix, the state transition probability difference D between any two hidden states i and j is calculated. trans (i, j), the calculation formula is:

[0106] ;

[0107] Where K is the total number of hidden states, k is the summation index variable, and a ika represents the probability of transitioning from hidden state i to hidden state k in the average state transition probability matrix. jk a represents the probability of transitioning from hidden state j to hidden state k in the average state transition probability matrix. ki a represents the probability of transitioning from hidden state k to hidden state i in the average state transition probability matrix. kj The probability of transitioning from hidden state k to hidden state j in the average state transition probability matrix is ​​represented by ||, where || represents the absolute value operation.

[0108] In a further embodiment, determining the optimal number of hidden states further includes:

[0109] The state redundancy index R(K) is obtained by calculating the minimum difference in state transition probabilities among all hidden state pairs. The formula is R(K) = min i≠j D trans (i, j);

[0110] A joint optimization objective function is constructed based on the state redundancy index R(K) and the Bayesian information criterion. The optimal number of hidden states K* is determined by minimizing the joint optimization objective function.

[0111] K*=argmin K {BIC(K)-λ×log(R(K)+ε)};

[0112] Where λ is the balance coefficient; ε is the smoothing constant; and BIC(K) is the Bayesian information criterion.

[0113] Alternatively, the Bayesian information criterion BIC(K) can be calculated based on the parameter estimation results of the Hidden Markov Model; a joint optimization criterion is constructed based on the state redundancy index R(K) and the Bayesian information criterion BIC(K); and the optimal number of hidden states K* is determined by minimizing the joint optimization objective function. The calculation formula is: K* = argmin K {BIC(K)-λ×log(R(K)+ε)}; where λ is a preset balance coefficient used to control the weights of model fit and state discrimination; ε is a preset smoothing constant used to prevent the logarithmic function from diverging.

[0114] In this embodiment, using R(K) alone can only evaluate redundancy, not goodness of fit; using BIC(K) alone can only evaluate the balance between fit and complexity, but cannot avoid physical redundancy. A penalty term -λ×log(R(K)+ε) is introduced through a joint optimization criterion. Here, R(K) is the difference between the most similar pair of states in the current model, representing the worst discriminative ability of the model. When there are highly redundant states in the model, R(K) will approach 0, causing -log(R(K)) to become very large, increasing the value of the objective function, thus penalizing that K value and causing it to be rejected.

[0115] Furthermore, the parameter λ plays a crucial role in adjusting the weights for statistical fitting and physical discrimination. In some preferred embodiments, the value of λ ranges from [0.5, 2.0]. When λ is large, the model tends to select a concise model with extremely high state discrimination; when λ is small, the model tolerates a certain degree of state similarity in exchange for higher fitting accuracy. The parameter ε is a very small positive number, preferably set to 1×10. -6 Its function is solely to ensure the stability of mathematical calculations and prevent the logarithmic function from becoming meaningless when R(K)=0. By traversing K∈{2,3,4,5,...,N} to calculate the objective function value, the K* that minimizes the function value is selected as the optimal number of hidden states. This ensures that the final output reservoir operation state sequence accurately reconstructs historical hydrological processes and possesses clear and independent physical scheduling meaning. Here, N is a natural number greater than 0.

[0116] This embodiment addresses the problems of traditional AIC or BIC criteria, which only consider statistical goodness of fit and are prone to causing redundant states in the model with repeated physical meanings. It proposes a joint optimization strategy that combines physical discrimination and statistical goodness of fit.

[0117] like Figure 5 As shown, according to one aspect of this application, before determining the optimal number of hidden states, an input variable selection step based on variable importance feedback is further included, specifically:

[0118] Extract the variable importance scores for each feature variable in the decision tree model;

[0119] The multidimensional hydrological variables are sorted according to their importance scores, and redundant variables with scores below a preset threshold are removed to obtain the filtered variable combination.

[0120] Reconstruct the time series based on the filtered variable combinations and return to the steps of constructing the Hidden Markov Model.

[0121] In other words, the time series is reconstructed based on the filtered variable combinations, and the steps of building a hidden Markov model are returned until the variable importance scores of all retained variables are not lower than the preset threshold.

[0122] In this embodiment, during the initial construction of the HMM, all collected hydrological variables such as Q may be included. in Q outVariables such as WL, Volume, downstream flow, and rainfall are used as inputs. However, not all variables play a crucial role in classifying the reservoir's operational status. Introducing irrelevant variables may increase computational complexity and introduce noise. It is preferable to utilize the built-in property of CART decision trees, namely variable importance, for selection. Variable importance is typically measured based on the total decrease in the Gini coefficient caused by the variable during tree splitting.

[0123] The results of the variable importance ranking under different input combinations are as follows: Figure 6 As shown. Specifically, taking the actual data of Reservoir A as an example, after preliminary modeling and CART training, the extracted variable importance ranking results may show: WL (water level) ≈ Volume (reservoir capacity) > Q out (Outbound flow) > Q in (Inflow); where outflow is the discharge. This indicates that in the scheduling mode identification of Reservoir A, water level and reservoir capacity are the decisive state variables, corresponding to rules such as water storage and flood control restrictions, while outflow is secondary, and inflow is the least important, because inflow is a natural process and not directly controlled, while the scheduling state reflects more the intention of human control. Based on this, if Q is found in If the importance score is significantly lower than other variables, such as being below a preset threshold of 5% or having an order-of-magnitude difference from other variables, then Q is determined. in This is a relatively redundant variable. In the next iteration, Q will be removed. in Only use [WL, Volume, Q] out Reconstruct the time series and train the Hidden Model (HMM). Through a feedback loop mechanism, it can automatically converge to the optimal combination of input variables. For example, for reservoir A, the final optimal model input combination might be Q. out And Volume, or Q out The combination of WL and Volume. It should be noted that the preset threshold can be determined based on the total number of feature variables and the distribution characteristics of variable importance scores.

[0124] This embodiment provides a post-model feature engineering method that uses the interpretability of decision trees to feed back into the construction of the hidden Markov model, thereby eliminating redundant noise variables and improving the robustness of state recognition.

[0125] In another embodiment of this application, homogeneous HMM remains an effective and computationally inefficient basic solution for simple reservoir scenarios with limited data, scarce computing resources, or indistinct seasonality.

[0126] Specifically, when constructing the Hidden Markov Model, it is assumed that the transition probabilities of the reservoir's operating states do not change with time, i.e., the state transition probability matrix A is a constant matrix. For any time t, the probability P(Si) of transitioning from state i to state j is given by...t= j|S t-1 =i)=a ij The parameters remain constant. The parameter estimation process employs the standard Baum-Welch algorithm, also known as the standard EM algorithm. In the expectation step, the forward-backward algorithm is used to calculate the statistics. At this point, the forward variable α... t The recursive formula for (j) is:

[0127] ;

[0128] In the maximization step, the transition probability parameters are updated directly by normalizing the frequency counts:

[0129] ;

[0130] Where, ξ t (i, j) represents the expected number of times that a character is in state i at time t and in state j at time t+1, γ t (i) represents the expected number of times the state is in state i at time t, and T is the total length of the time series.

[0131] Compared to non-homogeneous models, this embodiment reduces the number of model parameters to only K*(K-1) transition parameters and eliminates the need to estimate the regression coefficient β, resulting in faster training and more stable convergence. However, this comes at the cost of sacrificing the ability to explicitly model the dynamic switching mechanism between flood season and non-flood season. In practical applications, the method in this embodiment can be used as a benchmark model for comparison with non-homogeneous models, or as an alternative when covariate data is missing, such as when accurate inflow rates cannot be obtained.

[0132] In summary, a method for identifying reservoir operation status by coupling a hidden Markov model and a decision tree includes: acquiring multidimensional hydrological observation variables of the reservoir to construct a time series; constructing a non-homogeneous hidden Markov model, introducing physical factors reflecting the reservoir's hydrological situation as time covariates, establishing a dynamic mapping relationship between state transition probabilities and time covariates, and obtaining a hidden state sequence adapted to the seasonal characteristics of reservoir scheduling through parameter estimation and decoding; using the hidden state sequence as labels to train a decision tree model, extracting time-segmented reservoir operation discrimination rules; and constructing a model optimization index based on the difference in transition probabilities between states to determine the optimal number of hidden states.

[0133] According to one aspect of this application, a reservoir operation status identification device coupled with a hidden Markov model and a decision tree is provided, comprising:

[0134] Memory is used to store computer programs.

[0135] In this embodiment, the memory may include volatile memory and non-volatile memory, wherein the volatile memory may be random access memory (RAM); and the non-volatile memory may be read-only memory (ROM), flash memory, hard disk drive (HDD), or solid-state drive (SSD). The memory serves as the physical carrier of data and is used to store the operating system, device drivers, and core business logic code.

[0136] Specifically, the memory stores computer program instructions specifically designed to implement the reservoir operation status identification method that couples a hidden Markov model with a decision tree. These instructions contain the coding logic for each step in the aforementioned embodiments. Furthermore, the memory also temporarily stores the multidimensional hydrological variable time series generated during processing, the calculated time-varying transition probability matrix A(t), the hidden state sequence results, and the extracted discrimination rule set.

[0137] The processor is configured to execute a computer program stored in the memory, and when the computer program is executed, implement the reservoir operation status identification method coupled with a hidden Markov model and a decision tree as described in any of the above embodiments.

[0138] Specifically, the processor is the core of the device's computing and control system, and can be one or more central processing units (CPUs), digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices. The processor communicates with the memory via a bus, reading and decoding computer program instructions from the memory. When the processor executes a program, it performs the following operations: controlling the input interface to acquire reservoir operation observation data; and calling the mathematical calculation unit to perform complex matrix operations, including calculating the linear predictor η based on covariates. ij (t) Run the extended EM algorithm for parameter iterative optimization and execute the time-varying Viterbi decoding algorithm; use the logic operation unit to construct and split the CART decision tree; determine the optimal model structure based on the calculated state redundancy index R(K) and joint optimization criteria, and feed the identification results back to the scheduler through the output interface, such as the display screen or network port.

[0139] In some practical applications, the device can be manifested as a server deployed in a reservoir dispatch center, a high-performance workstation, or an edge computing terminal that integrates hydrological monitoring functions.

[0140] In one possible implementation, when constructing the non-homogeneous hidden Markov model, the time covariate X is... t Treating them as exogenous variables means assuming that the observed values ​​of the covariates are independent of the hidden state S at the current time. t generate.

[0141] In this embodiment, the time covariate Xt Includes inflow factor Q in (t), outflow factor Q out (t), water level factor WL(t), and reservoir capacity factor Volume(t). From a physical mechanism perspective, the assumption that these variables are exogenous variables is reasonable in the following ways: inflow factor Q in (t) is determined by the upstream rainfall and runoff process, and belongs to the uncontrollable natural inflow, which is not affected by the current reservoir operation status; the water level factor WL(t) and the storage capacity factor Volume(t) are the cumulative results of the water balance at the previous moment, and their values ​​are mainly determined by the storage and release process at time t-1 and before, which have become the established boundary conditions when observed at time t; the outflow factor Q out Although affected by scheduling decisions, at the daily modeling scale, the outflow for the current day is usually scheduled in the early morning of the current day or the evening of the previous day, which can be regarded as a known prerequisite for the state transition of the current day. Therefore, the covariate vector X at time t is used. t Calculating the state transition probability from time t-1 to time t is logically sound. This assumption is also consistent with the classic application of non-homogeneous hidden Markov models in speech recognition and bioinformatics, where covariates are input to the model as external driving signals, rather than being endogenous products of the hidden states.

[0142] In an exemplary embodiment, hydrological observation data from a typical day during the flood season at Reservoir A were substituted into a multinomial Logistic regression model for calculation. The resulting state transition probabilities show that under high water level and high inflow conditions, the probability of the reservoir transitioning from its current operational state is the highest; the probability of transitioning to a high-flow-rate flood discharge state is lower but significantly higher than during the dry season; and the probability of transitioning to a water storage state is extremely low. The results are consistent with the physical law that reservoirs maintain a stable operational mode and cease water storage during the flood season, verifying the effectiveness of covariate-driven transition probabilities.

[0143] According to one aspect of this application, when constructing a multinomial Logistic regression model, in addition to setting an identifiability constraint that the linear predictor of the reference state is zero, a regularization constraint can also be introduced as needed to prevent parameter overfitting.

[0144] Specifically, since the product of covariates and regression coefficients constitutes a linear predictor, multicollinearity among covariates, such as a high correlation between water level WL and reservoir capacity, can lead to unstable estimations of regression coefficients. Therefore, an L2 regularization penalty term is introduced during the parameter update process in the maximization step.

[0145] Specifically, the original weighted log-likelihood function is modified into a regularization objective function L. reg (β):

[0146] ;

[0147] Where L(β) is the original weighted log-likelihood function, μ is the regularization coefficient, and β ij，m Let ∑ be the regression coefficient of the m-th covariate when transitioning from state i to state j. i，j，m This represents the summation of regression coefficients across all non-reference states.

[0148] The regularization coefficient μ was determined through cross-validation. Specifically, the time series was divided into training and validation sets by year, and the value of μ that maximized the log-likelihood of the validation set was selected from the candidate value set {0.001, 0.01, 0.1, 1.0}. For the case data of Reservoir A, the optimal regularization coefficient determined through cross-validation was μ = 0.01.

[0149] Furthermore, physical constraints are set for the range of values ​​for the regression coefficients.

[0150] In this embodiment, based on the physical laws of reservoir scheduling, sign constraints or numerical boundaries can be set for some regression coefficients. For example, for the transition probability from an arbitrary state to a high-flow-rate flood discharge state, the inflow factor Q... in The regression coefficient should be positive, i.e., β. i3，1 A value greater than 0 reflects the physical logic that a larger inflow rate is more likely to trigger flood discharge scheduling; the regression coefficient of the water level factor WL for the transition to the water storage state should be negative, i.e., β. i1，3 A value less than 0 reflects the scheduling constraint that higher water levels make further water storage less suitable. In the maximization step of parameter estimation, boundary constraints can be implemented using constrained optimization algorithms to limit the regression coefficients to a physically reasonable range.

[0151] According to another aspect of this application, for non-homogeneous hidden Markov models, a time-averaging strategy is used to convert the time-varying state transition probability matrix sequence A(t). t=1 T It is transformed into a single representative transition probability matrix.

[0152] In this embodiment, due to the state transition probability α of the non-homogeneous hidden Markov model ij (t) changes dynamically with time covariates, making it unsuitable for directly calculating the difference in transition probabilities between states that measures the distinction between any two hidden states in their transition behaviors. To address this issue, a time-averaged transition probability matrix A* is defined, with its elements calculated using the following formula:

[0153] ;

[0154] Where T is the total length of the time series, a ij (t) represents the probability of transitioning from state i to state j at time t, where a* ij This represents the average transition probability.

[0155] The average transition probability reflects the overall tendency of the reservoir to switch between states throughout the entire observation period, eliminating the interference of extreme values ​​of covariates at a single moment on the calculation of the degree of difference, and obtaining a robust estimate of the dynamic relationship between states.

[0156] In another detailed embodiment, the transformation of time-varying transition probabilities to fixed values ​​and the difference in state transition probabilities D between states are demonstrated using a numerical case with K=4. trans A comprehensive demonstration of the complete calculation process.

[0157] Suppose a non-homogeneous hidden Markov model with K=4 is constructed based on daily data of Reservoir A from 2015 to 2023. The four hidden states represent State1 (dry season water storage), State2 (conventional water supply), State3 (flood season flood control), and State4 (ecological scheduling). After time averaging, the representative state transition probability matrix A* is as follows:

[0158] Line 1 (starting from State1): a* 11 =0.92, a* 12 =0.06, a* 13 =0.01, a* 14 =0.01.

[0159] Line 2 (starting from State2): a* 21 =0.08, a* 22 =0.85, a* 23 =0.04, a* 24 =0.03.

[0160] Line 3 (starting from State3): a* 31 =0.02, a* 32 =0.05, a* 33 =0.88, a* 34 =0.05.

[0161] Line 4 (starting from State4): a* 41 =0.03, a* 42 =0.07, a* 43 =0.05, a* 44 =0.85.

[0162] Now calculate the difference in state transition probabilities D between State1 and State2. trans (1, 2).

[0163] D trans (1, 2) = 1 / 2 × ∑ k=14 |a* 1k -a* 2k |+1 / 2×∑ k=1 4 |a* k1 -a* k2 |=1 / 2×1.68+1 / 2×1.70=1.69.

[0164] Similarly, calculate the dissimilarity of other state pairs:

[0165] D trans (1, 3) = 1 / 2 × (|0.92 - 0.02| + |0.06 - 0.05| + |0.01 - 0.88| + |0.01 - 0.05|) + 1 / 2 × (|0.92 - 0.01| + |0.08 - 0.04| + |0.02 - 0.88| + |0.03 - 0.05|) = 1.825.

[0166] D trans (1, 4) = 1 / 2 × (|0.92 - 0.03| + |0.06 - 0.07| + |0.01 - 0.05| + |0.01 - 0.85|) + 1 / 2 × (|0.92 - 0.01| + |0.08 - 0.03| + |0.02 - 0.05| + |0.03 - 0.85|) = 1.795.

[0167] D trans (2,3)=1 / 2×(|0.08-0.02|+|0.85-0.05|+|0.04-0.88|+|0.03-0.05|)+1 / 2×(|0.06-0.01|+|0.85-0.04|+|0.05-0.88|+|0.07-0.05|)=1.715.

[0168] D trans (2, 4) = 1 / 2 × (|0.08 - 0.03| + |0.85 - 0.07| + |0.04 - 0.05| + |0.03 - 0.85|) + 1 / 2 × (|0.06 - 0.01| + |0.85 - 0.03| + |0.05 - 0.05| + |0.07 - 0.85|) = 1.655.

[0169] D trans (3, 4) = 1 / 2 × (|0.02 - 0.03| + |0.05 - 0.07| + |0.88 - 0.05| + |0.05 - 0.85|) + 1 / 2 × (|0.01 - 0.01| + |0.04 - 0.03| + |0.88 - 0.05| + |0.05 - 0.85|) = 1.65.

[0170] Summarize the differences between all state pairs and calculate the state redundancy index R(K). Specifically, when K=4, there are 6 state pairs, and the differences between each pair are summarized as follows:

[0171] D trans (1, 2) = 1.69; D trans (1, 3) = 1.825; D trans (1, 4) = 1.795; D trans (2, 3) = 1.715; D trans (2, 4) = 1.655; D trans (3, 4) = 1.65.

[0172] The state redundancy index R(K) is equal to the minimum difference between all state pairs:

[0173] R(4)=min{1.69, 1.825, 1.795, 1.715, 1.655, 1.65}=1.65.

[0174] The minimum value corresponds to the difference between State3 and State4, indicating that in the K=4 model, the transition behaviors of flood control state and ecological scheduling state are most similar, and they are potentially redundant state pairs. The value of R(4)=1.65 is relatively high (the theoretical maximum value is 2), indicating that the state discrimination of the K=4 model is generally good.

[0175] In one embodiment of this application, the dimensional difference between BIC(K) and log(R(K)+ε) can be standardized by using the Min-Max normalization method.

[0176] In this embodiment, the value of the Bayesian information criterion BIC(K) is typically in the range of thousands to tens of thousands, while when R∈[0,2] and ε=10 -6 At this time, the value of log(R(K)+ε) ranges from approximately -13.8 to 0.69. To ensure that the two terms have comparable weight contributions in the objective function, they are normalized separately. Let the set of candidate hidden states be K∈{2, 3, 4, 5, 6}, calculate BIC(K) and R(K) for each K value, and then perform the following normalization:

[0177] BIC norm (K)=(BIC(K)-BIC min ) / (BIC max -BIC min );

[0178] R norm (K)=(R(K)-R min ) / (R max -R min).

[0179] Among them, BIC min and BIC max Let R be the minimum and maximum BIC values ​​among all candidate K values. min and R max These are the minimum and maximum values ​​of the state redundancy index R among all candidate K values, respectively, and BIC. norm (K) is the normalized Bayesian information content criterion, R norm (K) is the normalized state redundancy index. After normalization, BIC norm (K) and R norm (K) all fall within the interval [0, 1], eliminating the influence of dimensional differences.

[0180] In a further embodiment, an adaptive method based on data features is used to determine the optimal value of the balance coefficient λ.

[0181] In this embodiment, the balance coefficient λ is used to control the relative weights of statistical goodness of fit and state discrimination in the objective function. The optimal value of λ is related to the data characteristics of the specific reservoir and is determined using the following adaptive method: The candidate value set for λ is set as {0.5, 1.0, 1.5, 2.0, 2.5, 3.0}; for each candidate λ value, the corresponding optimal number of hidden states K*(λ) is determined according to the joint optimization objective function; for each model corresponding to K*(λ), the mean log-likelihood value on the validation set is calculated using the one-year leave-one-year cross-validation method; the λ value that maximizes the mean log-likelihood of the validation set is selected as the final balance coefficient. For the case of Reservoir A, the optimal balance coefficient determined through the adaptive process is λ=1.5. In some simplified application scenarios, λ=1.0 can also be directly used as the default value, indicating that statistical goodness of fit and state discrimination are equally important.

[0182] In another embodiment of this application, a method for identifying the operational status of a reservoir by coupling a Hidden Markov Model (HMM) and a Decision Tree (CART) model is provided. Based on conventional hydrological time-series data observed in the reservoir, the method specifically includes the following steps:

[0183] Step 1: Data collection and preprocessing.

[0184] Hydrological data and basic information about the target reservoir are collected. To ensure the accuracy and stability of model training, the data needs to be preprocessed, including outlier identification and removal, and missing value imputation. The processed data forms a continuous time series matrix, which serves as the data foundation for subsequent model training.

[0185] The hydrological data includes daily inflow and inflow Q. in Outbound flow Qout Key variable data closely related to the reservoir's operational status include water volume and reservoir level (WL). Basic reservoir information data includes the reservoir's catchment area, normal water level, dead water level, normal water level reservoir capacity, dead storage capacity, regulating capacity, installed capacity, and development tasks.

[0186] Missing data imputation, such as reservoir outflow or inflow, is primarily achieved by fitting the flow data from the nearest river hydrological station. Specifically, if a reservoir lacks a segment of outflow, and there are no significant tributaries flowing into the river between the nearest downstream hydrological station and the reservoir dam site, then the fitting relationship between the time series data of other outflows that are not missing and the corresponding flow data from the river hydrological station is analyzed to impute the missing outflow data. Furthermore, due to the characteristics of the Hidden Markov Model (HMM), if the time series data length is less than 100 days, it is best to discard that data; that is, the data length used in the model should be at least greater than 100 consecutive days.

[0187] Based on the preprocessed data, an inbound flow Q is generated. in Outbound flow Q out The time series matrix of multidimensional observation variables, including effective water storage volume and reservoir water level WL.

[0188] Step 2: Identify potential operating states based on Hidden Markov Models.

[0189] A Hidden Markov Model (HMM) is constructed based on the obtained reservoir time series matrix. The model structure includes: input variables Vars, number of hidden states K, response model, prior model, and transition model. The number of hidden states is set to K (K≥2), with each state representing a typical operating mode, such as storage period, flood discharge period, and ecological discharge. The probability density distribution type of hydrological variables is identified based on their distribution characteristics. A response model is constructed based on the distribution of observed variables. The prior model and transition model are constructed using the `transInit()` function of the `depmixS4` package. The HMM model is then constructed using the `makeDepmix()` function of the `depmixS4` package. During training, the model parameters, including the state transition probability matrix and the probability distribution parameters of observed variables, are iteratively updated using the `fit()` function of the `depmixS4` package until the log-likelihood convergence or the set maximum number of iterations is reached. If the model fails to converge, the initial parameters or the number of hidden states can be adjusted appropriately. Using the Viterbi decoding algorithm, the observation sequence at each time step is decoded to obtain the most probable hidden state label sequence S. tThis sequence not only characterizes the temporal changes in the reservoir's operational status, but also reflects the persistence and transition characteristics of the status.

[0190] For response model construction, it is necessary to first determine the distribution of the observed variables, and then construct the response model based on the distribution of the observed variables. Specifically, the hydrological data of the reservoir is cleaned to remove missing values. The distribution distribution is plotted by combining the distribution fitting enhancement package (fitdistrplus) in R language and the Akaike Information (AIC) criterion. The distribution type of the data is determined based on the best fit.

[0191] At daily or higher time resolutions, low flow rates occur frequently, while high flow rates (such as flood events) occur with low probability but large amplitude. Therefore, flow rate data distributions generally exhibit significant long-tail or skewed characteristics. Reservoir water levels and storage capacities are constrained by reservoir scheduling rules, resulting in high concentrations near the controlled water level and corresponding storage capacity ranges, generally showing obvious non-normal characteristics. The data distribution types in the depmixS4 package include normal, Gaussian, and Poisson distributions, which do not match the distribution characteristics of hydrological time series such as flow rate, water level, and storage capacities. For commonly used distribution types in hydrological time series such as log-normal, gamma, and Weiber distributions, since these are not defined in the depmixS4 package, custom data distribution types are used to construct log-normal, gamma, and Weiber distributions to match the distribution characteristics of hydrological observation variables and improve model accuracy. The custom data distribution types are subsequently coupled with the HMM model through code block calls. Each distribution parameter was initialized using the maximum likelihood method and the method of moments estimation.

[0192] Step 3: Construct a decision tree model and extract the discrimination rules for the running state.

[0193] The optimal hidden state label sequence S output by the Hidden Markov Model t The target classification label is merged with the original input variable to form a supervised learning sample set for training the decision tree model. A classification method is used to train the decision tree model, with the Gini index as the splitting criterion. Hyperparameters such as the maximum depth, minimum number of leaf node samples, minimum split information gain, and minimum number of split samples are set to prevent overfitting. K-fold cross-validation is used to evaluate the model's stability and generalization ability, ensuring stable classification performance. After training, all classification paths of the decision tree are extracted and converted into running state discrimination rules to interpret the hidden state labels.

[0194] Step 4: Model evaluation and optimization.

[0195] To improve the generalization ability and classification accuracy of the coupled model, the constructed HMM-CART coupled model was evaluated and optimized. Specifically, this involved: constructing model input combinations with different numbers of hidden states and different variables; modeling different model input combinations; evaluating model performance and robustness based on model evaluation metrics, and selecting the optimal model.

[0196] The model input combinations include different numbers of hidden states and different variables. Increasing the number of hidden states leads to higher model complexity, and its upper limit is related to the reservoir development tasks and basic functions. The more functions the reservoir has, the more hidden functions are recommended to have. Model evaluation metrics include AIC (Akaike Information Criterion), BIC (Bayesian Information Criterion), log-likelihood, and accuracy. The calculation formulas are as follows:

[0197] AIC = 2K-2 Loglike;

[0198] BIC = ln(n) × K⁻² Loglike;

[0199] ;

[0200] Where K is the number of hidden states in the HMM model, n is the total number of samples used to train the model, and Loglike is the model's log-likelihood value. Smaller AIC and BIC values, and larger Loglike values, indicate better model performance. Accuracy is the precision between the state sequence obtained from the decision tree model and the state sequence extracted by the HMM model. S HMM，i The hidden state label identified by the HMM model for the i-th sample; S CART，i Φ is the hidden state label identified by the CART model for the i-th sample; Φ() is an indicator function, which takes the value 1 when the hidden state labels obtained by the two models are consistent, and 0 otherwise. The smaller the AIC and BIC values ​​and the larger the Loglike value, the better the model performance. The larger the Accuracy value, the closer the classification results of the Hidden Markov Model and the Decision Tree Model are, and the more reliable the model is.

[0201] For model input combinations (model input groups with different hidden state numbers × different variables), simulating all combinations exhaustively would consume unnecessary time and computational resources. In the HMM-CART model coupling process, a variable importance feedback mechanism based on the CART model is introduced to filter the input combinations of the HMM model. The initial identification results yield a ranking of variable importance, which is used to eliminate redundant variables and reconstruct the model's input combinations, thereby improving simulation and optimization efficiency.

[0202] Step 5: Model application and refined identification of reservoir operation status.

[0203] Based on the trained HMM-CART coupled model, inputting real-time or historical hydrological data time series from the reservoir and running the HMM-CART coupled model will output the reservoir operation status label corresponding to each time step. The final output includes the time series data of the operation status label and the rule paths corresponding to different operation statuses.

[0204] This invention employs a covariate-driven non-homogeneous modeling approach. By introducing physical factors such as water level, inflow, outflow, and reservoir capacity as time covariates, and utilizing multinomial logistic regression to construct a dynamic mapping relationship between these factors and state transition probabilities, the state transition probability matrix is ​​no longer constant but evolves in real time with hydrological conditions. This breaks the homogeneity constraint, enabling the model to accurately capture the dynamic switching logic under different operating conditions, such as flood control during the flood season and water storage during the non-flood season. This enhances the model's adaptability and physical interpretability to complex scheduling environments. Furthermore, a joint optimization approach based on physical difference is adopted. By defining the difference in transition probabilities between states, the physical distance between different states in their dynamic behavior of transitioning in and out is quantified. Based on this, a joint optimization criterion including a redundancy penalty term is constructed. While ensuring statistical goodness of fit, this criterion forcibly penalizes states with highly similar physical behaviors, effectively eliminating redundant states and ensuring that the finally identified operating states have clear and independent physical scheduling meanings.

[0205] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A method for identifying the operational status of a reservoir by coupling a hidden Markov model and a decision tree, characterized in that, include: Obtain reservoir operation observation data of the target reservoir and construct a time series containing multidimensional hydrological variables; A hidden Markov model is constructed, and the parameters of the hidden Markov model are estimated based on the time series. The time series is then decoded using the model with the estimated parameters to obtain the hidden state sequence that represents the reservoir operation mode. Using the hidden state sequence as the label and the time series as the feature, a decision tree model is constructed and trained to extract the reservoir operation discrimination rules corresponding to different hidden states. The optimal number of hidden states is determined based on the preset model optimization index, and the hidden state sequence corresponding to the optimal number of hidden states and the reservoir operation discrimination rule are output. The Hidden Markov Model (HMM) is a non-homogeneous HMM. The HMM contains a time-varying state transition probability matrix, and the state transition probabilities in the time-varying state transition probability matrix are functions of the time covariates. The time covariates should include at least the inflow factor Q, which reflects the characteristics of watershed inflow. in (t), the outflow factor Q, which reflects the intensity of reservoir discharge. out (t), water level factor WL(t) reflecting the current operating potential energy of the reservoir, and reservoir capacity factor Volume(t) reflecting the occupancy status of the reservoir; By establishing a dynamic mapping relationship between state transition probabilities and time covariates, the process of identifying hidden state sequences can be adapted to the seasonal and periodic characteristics of reservoir scheduling. A dynamic mapping relationship is constructed using a multinomial logistic regression model. The specific calculation process includes: Define the time covariate vector X at time t. t For Q in (t), Q out Combinations of (t), WL(t), and Volume(t); Calculate the linear predictor η for the transition from hidden state i to hidden state j at time t. ij (t), the calculation formula is: ； Where, β ij，0 For the baseline intercept term, β ij，1 To β ij，4 These are the regression coefficients for inflow factor, outflow factor, water level factor, and reservoir capacity factor, respectively. Based on linear predictor η ij (t), calculate the state transition probability P(S) from hidden state i to hidden state j at time t. t =j|S t-1 =i,X t The calculation formula is: ； Among them, S t S represents the hidden state at time t. t-1 Let K represent the hidden state at time t-1, and K be the total number of hidden states. x It is an exponential function, and k is the summation index variable.

2. The method according to claim 1, characterized in that, Constructing a Hidden Markov Model, including constructing a response model, specifically involves: Identify the probability density distribution characteristics of multidimensional hydrological variables, and construct response models using non-normal distribution functions for variables with extreme value distributions; Among them, the non-normal distribution function includes the Weibull distribution function and the log-normal distribution function. The Weibull distribution function is used to describe the probability distribution of the inflow variable, and the log-normal distribution function is used to describe the probability distribution of the outflow, reservoir water level and reservoir storage variables.

3. The method according to claim 1, characterized in that, Extract the reservoir operation discrimination rules corresponding to different hidden states, specifically including: The state labels in the hidden state sequence are merged with the corresponding multidimensional hydrological variables and time covariates to construct a supervised learning sample set. Based on the pre-defined reservoir scheduling cycle characteristics, the supervised learning sample set is divided into a flood season sample set and a non-flood season sample set; Pre-configured classification and regression tree models were trained for flood season and non-flood season sample sets respectively. The Gini index was used as the splitting criterion to extract the reservoir operation discrimination rules for different time periods.

4. The method according to claim 1, characterized in that, The optimal number of hidden states is determined based on preset model optimization indices, including calculating the difference in transition probabilities between states: Calculate the difference in state transition probabilities D between any two hidden states i and j. trans (i, j), the calculation formula is: ； Where K is the total number of hidden states, k is the summation index variable, and a ik Let a represent the probability of transitioning from hidden state i to hidden state k. jk Let a represent the probability of transitioning from hidden state j to hidden state k. ki Let a represent the probability of transitioning from hidden state k to hidden state i. kj Let | represent the probability of transitioning from hidden state k to hidden state j, and | represent the absolute value operation.

5. The method according to claim 4, characterized in that, Determining the optimal number of hidden states also includes: The state redundancy index R(K) is obtained by calculating the minimum difference in state transition probabilities among all hidden state pairs. The formula is R(K) = min i≠j D trans (i, j); A joint optimization objective function is constructed based on the state redundancy index R(K) and the Bayesian information criterion. The optimal number of hidden states K* is determined by minimizing the joint optimization objective function. K*=argmin K {BIC(K)-λ×log(R(K)+ε)}; Where λ is the balance coefficient; ε is the smoothing constant; and BIC(K) is the Bayesian information criterion.

6. The method according to claim 1, characterized in that, Before determining the optimal number of hidden states, an input variable selection step based on variable importance feedback is also performed, specifically: Extract the variable importance scores for each feature variable in the decision tree model; The multidimensional hydrological variables are sorted according to their importance scores, and redundant variables with scores below a preset threshold are removed to obtain the filtered variable combination. Reconstruct the time series based on the filtered variable combinations and return to the steps of constructing the Hidden Markov Model.

7. The method according to claim 1, characterized in that, Parameter estimation for the Hidden Markov Model involves using the Extended Expectation-Maximization (EMF) algorithm to solve for the regression coefficients of the multinomial Logistic Regression model, including: In the expectation step, based on the time series, the posterior probability of the hidden state is calculated using the forward-backward algorithm, where the state transition probability in the calculation process adopts the time-varying value obtained by the linear predictor at the current time. In the maximization step, a weighted log-likelihood function is constructed based on the posterior probability, and the regression coefficients are iteratively updated using the iterative weighted least squares method or the quasi-Newton method until the model converges.

8. A reservoir operation status identification device coupled with a hidden Markov model and a decision tree, characterized in that, include: Memory, used to store computer programs; A processor is configured to execute a computer program stored in the memory, wherein when the computer program is executed, it implements the reservoir operation status identification method coupled with a hidden Markov model and a decision tree as described in any one of claims 1 to 7.

Images