Lagrangian Relaxation Based on Improved Proximal Bundle Method for Short-Term Hydrothermal Scheduling

: Short-term hydrothermal scheduling (STHS) can improve water use efﬁciency, reduce carbon emissions, and increase economic beneﬁts by optimizing the commitment and dispatch of hydro and thermal generating units together. However, limited by the large system scale and complex hydraulic and electrical constraints, STHS poses great challenges in modeling for operators. This paper presents an improved proximal bundle method (IPBM) within the framework of Lagrangian relaxation for STHS, which incorporates the expert system (ES) technique into the proximal bundle method (PBM). In IPBM, initial values of Lagrange multipliers are ﬁrstly determined using the linear combination of optimal solutions in the ES. Then, each time PBM declares a null step in the iterations, the solution space is inferred from the ES, and an orthogonal design is performed in the solution space to derive new updates of the Lagrange multipliers. A case study in a large-scale hydrothermal system in China is implemented to demonstrate the effectiveness of the proposed method. Results in different cases indicate that IPBM is superior to standard PBM in global search ability and computational efﬁciency, providing an alternative for STHS.


Introduction
According to the International Energy Agency, thermal power and hydropower are basic sources of electricity production in many countries [1]. Thus, short-term hydrothermal scheduling (STHS) is necessary in power system operations. The significance of STHS is to improve water use efficiency, reduce carbon emissions, and increase economic benefits by optimizing the commitment and generation level of hydro and thermal generating units together [2]. However, limited by complex hydraulic and electrical constraints, the nature of STHS is a large-scale nonconvex, nonlinear problem with integer variables, posing great challenges in modeling for operators [3].
Many approaches have been developed for STHS, such as mixed-integer linear programming (MILP) [4], nonlinear programming (NLP) [5], dynamic programming (DP) [6], and genetic algorithm (GA) [7]. Although these approaches have achieved success in practice, challenges still exist when dealing with large-scale systems. MILP depends on linearization strategies to accurately represent the system behavior. For NLP, a convex approximation is required to improve the global search ability, and global optimum is not guaranteed. DP suffers from the "curse of dimensionality". GA falls into premature convergence easily. For large-scale systems, common approaches are based on decomposition techniques, such as the Lagrangian relaxation (LR) [8], Benders decomposition (BD) [9], and aggregation-disaggregation [10]. BD decomposes the primal problem into a master problem and some subproblems for dimension reduction. However, BD is sensitive to integer variables, of which the computational efficiency becomes very low for Although these plants belong to different generation companies, they are operated by a central dispatching center in a centralized manner. Moreover, according to local government policy, the operation objective of the hydrothermal system is to minimize coal consumption to achieve energy-saving and emission reduction. For such a large system, a small improvement in operations will translate into huge benefits, making STHS an important problem.

Objective Function
The objective function minimizes the coal consumption of the hydrothermal system over the whole scheduling horizon:  Although these plants belong to different generation companies, they are operated by a central dispatching center in a centralized manner. Moreover, according to local government policy, the operation objective of the hydrothermal system is to minimize coal consumption to achieve energy-saving and emission reduction. For such a large system, a small improvement in operations will translate into huge benefits, making STHS an important problem.

Optimization Model 2.2.1. Objective Function
The objective function minimizes the coal consumption of the hydrothermal system over the whole scheduling horizon: min p i,j,t ,p k,m,t ,u i,j,t ,u k,m,t where t, i, j, k, and m = indices of time period, thermal plant, thermal unit, hydro plant, and hydro unit, respectively; T = number of time periods; N TP = number of thermal plants; N i TU = number of thermal units in plant i; p i,j,t = power output of thermal unit j in plant i at period t (MW); p k,m,t = power output of hydro unit m in plant k at period t (MW); u i,j,t = binary variable indicating whether or not thermal unit j in plant i is operating at period t; u k,m,t = binary variable indicating whether or not hydro unit m in plant k is operating at period t; a i,j (t/MW2h), b i,j (t/MWh), and c i,j (t/h) = fuel consumption coefficients of thermal unit j in plant i, respectively; ∆t = time conversion variable (h).
It is important to point out that because hydro units and thermal units are closely coupled by electrical constraints (Equations (2) and (3)), the schedules of both hydro units and thermal units will affect the objective function. Thus, the decision variables in Equation (1) are p i,j,t , p k,m,t , u i,j,t , and u k,m,t .

Constraints
In STHS, constraints that should be satisfied are grouped into system-wide electrical constraints, thermal power constraints, and hydropower constraints, which are formulated as follows.
1. System-wide electrical constraints where N HP = number of hydro plants; N k HU = number of hydro units in plant k; D t = system load demand at period t (MW); p max i,j and p max k,m = the maximum output limits of thermal unit and hydro unit (MW), respectively; α = the constant representing the percentage of spinning reserve in load demand.
Equation (2) sets the power balance constraints. Equation (3) sets the spinning reserve constraints.

Thermal power constraints
where p min i,j = the minimum output limit of unit j in plant i (MW); y i,j,t = binary variable indicating whether or not the unit is started up at period t; z i,j,t = binary variable indicating whether or not the unit is shut down at period t; SU i,j and SD i,j = start-up and shut-down ramp rate limits of unit j in plant i (MW/h), respectively; RU i,j and RD i,j = ramp-up and ramp-down limits of unit j in plant i (MW/h), respectively; TA i,j and TB i,j = the minimum online and offline time periods of unit j in plant i (h), respectively. Equation (4) sets the lower and upper bounds of each unit on power output. Equations (5)-(7) set the power ramping constraints of each unit, including the start-up ramp rate limit, shut-down ramp rate limit, and normal ramp-up and ramp-down limits. Equations (8) and (9) denote the logical status of unit commitment. Equations (10) and (11) set the minimum online and offline time periods of each unit, respectively.

Hydropower constraints
where Q k,t , B k,t , and R k,t = total inflow, natural inflow, and total release of plant k at time t (m 3 /s), respectively; IU k = the direct upstream plants set of plant k; kh = the direct upstream plant index of plant k; τ kh,k = water travel time from plant kh to plant k (h); S k,t = storage of plant k at time t (m 3 ); R k,t and R k,t = total power release and spillage of plant k at time t (m 3 /s), respectively; R k,m,t = power release of unit m in plant = the minimum and maximum storages of plant k (m 3 ), respectively; Z k,t = forebay water level of plant k at time t (m); f k zs (·) = the function between storage and water level of plant k; ZT k,t = tailrace water level of plant k at time t (m); f k zr (·) = the function between tailrace water level and total release of plant k; H k,t = net water head of plant k at time t (m); f k hr (·) = the function between head loss and power release of plant k; S k,beg and S k,end = the initial and terminal storages of plant k (m 3 ), respectively; f k,m pqr (·) = power output as a function of unit power release and net head of unit m in plant k; p min k,m and p max k,m = the minimum and maximum power outputs of unit m in plant k (MW), respectively; ps max k,m,vz (H k,t ) and ps min k,m,vz (H k,t ) = the upper and lower bounds of vibration zone vz of unit m in plant k (MW), respectively. Both ps max k,m,vz (H k,t ) and ps min k,m,vz (H k,t ) are related to the net water head H k,t . Equations (12) and (13) set the water balance constraints. Equation (14) defines the total water release. Equations (15) and (16) set the discharge limit and storage limit, respectively. Equations (17) and (18) are forebay water level and storage function, and tailrace water level and release function, respectively, both of which are nonlinear and nonconvex. Equations (19) and (20) define the expression of the water head, where Equation (20) is nonlinear and convex. Equations (21) and (22) restrict the initial and terminal storages of each plant, respectively. Equation (23) is the hydropower production function, which is composed of a family of curves illustrating the nonlinear relationship between power output, water head, and power release. Equation (24) denotes the power output limit. Equation (25) sets the vibration zone limit.

Lagrangian Relaxation Framework
Within the LR framework, the first step is to construct a specific dual problem by relaxing linking constraints. Then, the next step consists of solving the dual problem. The obtained dual solution acts as a lower bound to the primal problem, and heuristics are adopted to convert the infeasible dual solution to a feasible solution. Through gradually narrowing the gap between dual solution and primal solution, a near-optimal solution will be finally found in multiple iterations.
In the context of STHS, Equations (1)- (25) are denominated the primal problem. The hydropower part and thermal power part are coupled by Equations (2) and (3). Therefore, the following dual function is constructed by relaxing Equations (2) and (3) via Lagrange multipliers λ D t and λ R t [28,29]: Subject to: Equations (4)- (25). After regrouping relevant terms, independent thermal and hydro subproblems can be split from Equation (26). Detailed forms of these subproblems can be found in [15]. To solve the thermal and hydro subproblems, MILP is adopted [30]. The reasons for using MILP are as follows: (1) the nonlinearity in hydro and thermal subproblems can be eliminated with piecewise linear techniques; (2) MILP is efficient in solving small-scale problems; and (3) MILP is easily programmed based on general mathematical solvers.
Moreover, the heuristic technique in [31] is adopted to transform the dual solution into a feasible solution, which fixes the hydropower solution first and then adjusts the thermal unit schedule using a priority-list approach. Finally, the iteration of LR stops if any of Equations (27) and (28) are satisfied: where n = the iteration index; rdg n , F * n and LD * n = the relative duality gap, the primal value (t), and the Lagrangian dual value (t) at iteration n, respectively; RDG = the supplied minimum relative duality gap; λ n D t and λ n R t = the values of λ D t and λ R t at iteration n, respectively; ∆λ = the supplied minimum multiplier variation.

Standard Proximal Bundle Method
In the dual function (Equation (26)), subgradients with respect to λ D t and λ R t are Let n be the current iteration index. The notation is simplified by denoting λ = λ D 1 , · · · , λ D T , λ R 1 , · · · , λ R T and g = g D 1 , · · · , g D T , g R 1 , · · · , g R T . PBM generates two , where nk = past iteration index; λ nk and g nk = the multiplier vector and subgradient vector at iteration nk, respectively; LD λ nk = the dual value associated with λ nk ; λ nk = the best multiplier vector so far at iteration nk, which provides the maximum dual value LD λ nk . The sequence LD λ nk , g nk nk<n is the so-called "bundle", while the sequence λ nk nk<n denotes the stability centers in the bundle [32].
Having the bundle, a cutting-planes model is formulated at iteration n LD n (λ) = min Then, the next iterate λ n+1 is obtained by solving the quadratic programming [33] where µ n = the penalty parameter, controlling the distance from λ n+1 to λ n . With the obtained λ n+1 , the ascent condition is checked: where ε = the parameter defining the minimum increase in dual value; δ n+1 = LD λ n+1 − LD λ n measuring the increase predicted by the cutting-planes model.
Although PBM is known for its stability and accuracy, there are still some drawbacks: (1) The time cost of solving quadratic programming. A quadratic programming problem in the form of Equation (31) needs to be solved at each iteration to generate a new iterate. The scale of the quadratic programming increases with the size of the bundle, which could be quite time-consuming in later iterations. (2) As the iteration proceeds, the ascent condition is hard to satisfy. As a result, many null steps are declared, causing PBM to easily fall into local optimum. (3) Initial values of Lagrange multipliers have an influence on the computational efficiency of PBM, which has not been fully discussed in the literature [33].

Improved Proximal Bundle Method
To overcome the difficulties of PBM, an IPBM incorporating the ES technique with PBM is proposed in this section. The ES component in IPBM mines information from historical operational records of hydrothermal units to guide the solution of a current problem. The ES consists of two parts: a knowledge base and inference engine, where the knowledge base is a repository of facts storing the knowledge about the STHS problem domain, and the inference engine provides reasoning about the information in the knowledge base to find a solution [25].
In the solution of Lagrangian dual, the main task is to determine the optimal values of Lagrange multipliers to reduce the relative duality gap (RDG). Thus, the trajectories of RDG and multiplier values in the iteration are regarded as a kind of knowledge. By extracting knowledge expressions and building the knowledge base, decision support can be provided for STHS from the inference engine. Specifically, initial values of multipliers can be obtained by using a linear combination of optimal multiplier values in the knowledge base, providing a good lower bound to the primal problem. Moreover, when a null step is claimed in IPBM, the solution space of STHS can be narrowed based on the knowledge base. Then, OD is carried out in the solution space to derive new updates of Lagrange multipliers, which can not only avoid solving quadratic programming, but also increase the likelihood of reaching a global optimum. Hence, the difficulties that standard PBM encountered are alleviated in IPBM.
Details about the ES technique are presented in the following.

Knowledge Base
To build the knowledge base, knowledge expressions are extracted from PBM iterations of historical scenarios. Since different generation scenarios correspond to different iteration processes, scenarios are firstly quantified by eigenvectors. Then, considering the large number of historical generation scenarios, cluster analysis is necessary to avoid the "scenario explosion" and the noise interference to the reasoning stage. Finally, knowledge expressions are extracted from the PBM iterations of representative scenarios and are stored in the database. Detailed steps of building the knowledge base are described as follows.
(1) Representation of generation scenario The generation scenario is characterized by operational constraints. Through interviewing experienced operators and analyzing the form of Equations (2)-(25), constraints and their associated factors are summarized in Table 1. In Table 1, constraints related to unit type and reservoir features usually do not change with scenario, or change slightly, while the rest of the constraints are dynamic. Therefore, the generation scenario is characterized by the load demand curve, natural inflows, and initial and terminal water levels of hydropower plants.

Constraint Associated Factor
(2), (3) Load demand curve (4)- (11) Thermal unit type (12) Natural inflows (13)- (20) Reservoir features (21), (22) Initial and terminal water levels (23)- (25) Hydro unit type To quantitatively represent the generation scenario, the following eigenvector η is formulated: where η 1 = total energy demand (MWh); η 2 = peak load (MW); η 3 = peak-valley difference in the load demand curve; η 4 = the vector of basin natural inflows; η 5 = the vector of basin storage energy; hb = hydro basin index; B hb = total natural inflows of hydro basin hb over the whole planning horizon (m 3 /s); HB = the number of hydro basins; PE hb = average storage energy of hydro basin hb over the whole planning horizon (MWh). In Equation (33), elements η 1 , η 2 , and η 3 depict the characteristics of the load demand curve; η 4 denotes the natural inflows of basins; η 5 means the average storage energy, which is a function of the initial and terminal water levels of hydropower plants [34].
(2) Cluster analysis of historical generation scenarios To obtain representative scenarios from massive historical generation scenarios, cluster analysis is necessary. First, the eigenvectors are normalized to eliminate the influence of dimension. Then, k-means clustering and silhouette coefficient are adopted to achieve scenario clustering [35]. By maximizing the silhouette coefficient of the clustering results, an optimal clustering result can be obtained. In each cluster, scenarios are sorted in ascending order according to the Euclidean distances from the centroid, and the scenario that has the minimum distance is selected as the representative scenario of the cluster.
(3) Knowledge expression extraction For each representative scenario, standard PBM is used to solve the Lagrangian dual. Knowledge expressions are extracted from the PBM iterations in the form of Equation (36): where nk and nk max = the iteration index and maximum iteration index in the process, respectively; λ nk and rdg nk = the multiplier vector and relative duality gap at iteration nk, respectively. Finally, all knowledge expressions and scenario eigenvectors are saved in the database to form the knowledge base, as illustrated in Figure 2.
where 1 η = total energy demand (MWh); 2 η = peak load (MW); 3 η = peak-valley difference in the load demand curve; 4 η = the vector of basin natural inflows; 5 η = the vector of basin storage energy; hb = hydro basin index; hb B = total natural inflows of hydro basin hb over the whole planning horizon (m 3 /s); HB = the number of hydro basins; hb PE = average storage energy of hydro basin hb over the whole planning horizon (MWh). In Equation (33), elements 1 η , 2 η , and 3 η depict the characteristics of the load demand curve; 4 η denotes the natural inflows of basins; 5 η means the average storage energy, which is a function of the initial and terminal water levels of hydropower plants [34].
(2) Cluster analysis of historical generation scenarios To obtain representative scenarios from massive historical generation scenarios, cluster analysis is necessary. First, the eigenvectors are normalized to eliminate the influence of dimension. Then, k-means clustering and silhouette coefficient are adopted to achieve scenario clustering [35]. By maximizing the silhouette coefficient of the clustering results, an optimal clustering result can be obtained. In each cluster, scenarios are sorted in ascending order according to the Euclidean distances from the centroid, and the scenario that has the minimum distance is selected as the representative scenario of the cluster.

Inference Engine
The inference engine works on the basis of the knowledge base. Similar generation scenarios mean that they have similar problem domains. Therefore, assuming that the

Inference Engine
The inference engine works on the basis of the knowledge base. Similar generation scenarios mean that they have similar problem domains. Therefore, assuming that the PBM iterations of similar scenarios are alike, initial values of Lagrange multipliers, and new updates of Lagrange multipliers when null steps occur, can be reasoned by the inference engine.
(1) Inferring the initial values of Lagrange multipliers The similarity between two scenarios can be measured by the Euclidean distance between their eigenvectors: where η * = eigenvector of scenario to be solved; η sn = scenario sn in the knowledge base.
Initial values of multipliers can be obtained by using a linear combination of optimal multiplier values of different scenarios. The weighting coefficient is inversely proportional to the Euclidean distance and is determined as follows: where SN = the set of scenarios in the knowledge base. Hence, the initial values of multipliers are obtained (Figure 3a): where λ sn is the optimal multiplier vector of scenario sn.
(1) Inferring the initial values of Lagrange multipliers The similarity between two scenarios can be measured by the Euclidean distance between their eigenvectors: where * η = eigenvector of scenario to be solved; sn η = scenario sn in the knowledge base. Initial values of multipliers can be obtained by using a linear combination of optimal multiplier values of different scenarios. The weighting coefficient is inversely proportional to the Euclidean distance and is determined as follows: where SN = the set of scenarios in the knowledge base. Hence, the initial values of multipliers are obtained (Figure 3a): where sn λ is the optimal multiplier vector of scenario sn . (2) Inferring new updates of Lagrange multipliers when null steps occur In the iteration process of PBM, when Lagrange multipliers fail to improve the value of objective function sufficiently, a null step will be declared. Assume that a null step occurs at iteration nu ; the relative duality gap is nu rdg ; and the most similar scenario in the knowledge base to the scenario to be solved is scenario se . The RDG measures the optimality of a dual solution. If two iterations have a close RDG, this means that they are at a (2) Inferring new updates of Lagrange multipliers when null steps occur In the iteration process of PBM, when Lagrange multipliers fail to improve the value of objective function sufficiently, a null step will be declared. Assume that a null step occurs at iteration nu; the relative duality gap is rdg nu ; and the most similar scenario in the knowledge base to the scenario to be solved is scenario se. The RDG measures the optimality of a dual solution. If two iterations have a close RDG, this means that they are at a similar stage in the whole iteration process. Therefore, rdg nu is used to match the knowledge expression in the knowledge base, as follows: where nk * = the iteration index in scenario se that has the closest RDG to rdg nu ; NK = the maximum iteration index in the iterations of scenario se.
With the obtained iteration index nk * , the range of multiplier values is derived: where λ max l and λ min l = the upper and lower bounds of the lth element in vector λ, respectively; λ nk l = the value of the lth element in vector λ at iteration nk. The obtained range is also the solution space of Lagrangian dual. Figure 3b depicts the procedure to infer new updates of Lagrange multipliers when null steps occur.
Having the solution space, a complete enumeration of state combinations of multiplier values can guarantee an optimum. However, the enumeration suffers from the "curse of dimensionality". Supposing that each of the multipliers is discretized into A states, there will be A 2T state combinations to deal with, which is hard to deal with. Thus, the concept of orthogonal design (OD) is introduced. The OD is an experimental design method for multiple-factor experiments, which samples a representative subset from complete combinations based on the orthogonal array [36]. For an experiment of X factors and Y levels per factor, there exists an orthogonal array shown as L C Y X = a i,j C×X , where L denotes the OD, Y X is the number of complete combinations, and C is the number of combinations to be tested in the OD. L C Y X is a matrix with C rows and X columns, and each row represents a combination of levels. Due to the orthogonality, the array L C Y X ensures that the sample combinations are scattered uniformly in the state space. Figure 4 depicts the difference between complete and orthogonal combinations of an experiment that has three factors with three levels per factor. It can be observed that the orthogonal combinations are representative and uniformly distributed in the state space. Moreover, the scale of orthogonal combinations is much smaller than complete combinations. similar stage in the whole iteration process. Therefore, nu rdg is used to match the knowledge expression in the knowledge base, as follows: The obtained range is also the solution space of Lagrangian dual. Figure 3b depicts the procedure to infer new updates of Lagrange multipliers when null steps occur.
Having the solution space, a complete enumeration of state combinations of multiplier values can guarantee an optimum. However, the enumeration suffers from the "curse of dimensionality". Supposing that each of the multipliers is discretized into A states, there will be 2T A state combinations to deal with, which is hard to deal with. Thus, the concept of orthogonal design (OD) is introduced. The OD is an experimental design method for multiple-factor experiments, which samples a representative subset from complete combinations based on the orthogonal array [36]. For an experiment of X factors and Y levels per factor, there exists an orthogonal array shown as , where L denotes the OD, X Y is the number of complete combinations, and C is the number of combinations to be tested in the OD.
( ) X C L Y is a matrix with C rows and X columns, and each row represents a combination of levels. Due to the orthogonality, the array ensures that the sample combinations are scattered uniformly in the state space. Figure 4 depicts the difference between complete and orthogonal combinations of an experiment that has three factors with three levels per factor. It can be observed that the orthogonal combinations are representative and uniformly distributed in the state space. Moreover, the scale of orthogonal combinations is much smaller than complete combinations. Additionally, before implementing the OD, constraint violation needs to be checked. Additionally, before implementing the OD, constraint violation needs to be checked.
where κ = the minimum power deviation (MW). If Equation (42) holds, it means that the power balance constraint at period t is broken. If Equation (43) holds, it means that the spinning reserve constraint at period t is broken.
Lagrange multipliers corresponding to broken constraints can be regarded as experiment factors. By scattering multipliers corresponding to broken constraints into several levels, a multiple-factor experiment is formed. Accordingly, an OD is executed to yield sample combinations. Keeping multipliers corresponding to unbroken constraints fixed, new updates of multipliers are obtained by combining sample combinations of broken multipliers with fixed unbroken multipliers. In this way, the computational burden of solving quadratic programming is avoided. More importantly, the OD behaves like a "mutation operator" to increase the likelihood of reaching global optimum. Then, after solving hydro and thermal subproblems with the obtained multiplier updates, arrange the dual values in ascending order, and update the bundle and stability centers sequences. Continue to solve the dual problem by PBM until a next null step is declared or the iteration stops.

Procedures of the Proposed IPBM
Procedures of the proposed IPBM within the LR framework are summarized as follows: Step 1: Build the knowledge base. First, represent historical generation scenarios by eigenvectors. Then, implement cluster analysis with k-means clustering and silhouette coefficient to obtain representative scenarios. Finally, extract knowledge expressions from PBM iterations of representative scenarios and save them in the database.
Step 2: Relax linking constraints to form the Lagrangian dual by Equation (26).
Step 3: Solve the dual problem. Solve hydro and thermal subproblems by MILP with the given multipliers. Then, set n = n + 1. If n is equal to two, go to Step 4. (d) If the ascent condition (Equation (32)) holds, a serious step is declared and go to Step 4. Otherwise, a null step is declared and update the multipliers by inference engine. (e) Solve the hydro and thermal subproblems by MILP. Then, according to the results, arrange the dual values in ascending order, and update the bundle and stability center sequences.
Step 4: Primary recovery. Generate a feasible solution by the heuristic used.

Parameter Settings and Performance Metrics of IPBM
The whole planning horizon of STHS is 24 h, with 1 h for each period. To build the knowledge base, historical operational data from 2013 to 2019 are used. There are 2380 generation scenarios at the daily scale after data correction, which are then clustered into 33 representative scenarios. Knowledge expressions are extracted from PBM iterations of these 33 representative scenarios. The knowledge base is based on MySQL Database 5.7 and structured query language. The proposed IPBM is encoded in Java language and implemented on a PC-Intel@2.60GHz. The MILP subproblems are solved with Gurobi 8.1.1 optimization solver. The OD is implemented using SPSS Statistics. Parameters associated with the computational simulation are set as follows: RDG = 0.5%, ∆λ = 10 −4 , ε = 10 −2 , and κ = 10 −2 .
To evaluate the performance of IPBM, three metrics are adopted: (1) dual value, representing the global search ability of IPBM; (2) primal value, representing the objective optimality after repairing the infeasible dual solution by heuristics; (3) computational time, representing the computational efficiency of IPBM. Among these metrics, the dual value is a positively oriented metric-the larger the better-while the remaining two metrics are negatively oriented metrics-the lower the better.
Based on the above settings, two improvements in IPBM, inferring initial values of Lagrange multipliers and inferring new updates of multipliers when null steps occur, are firstly tested. Then, the robustness of IPBM is demonstrated in comparison with standard PBM in 12 different cases.

Performance Testing of IPBM
A typical case in the spring of 2019 is selected as the generation scenario to manifest the performance of IPBM. In this scenario, water levels of hydropower plants are at the descending period before flood season. Both the load demand and water inflows are at a low level. The energy demand and peak load demand are 2.95 × 10 5 MWh and 1.51 × 10 4 MW, respectively. The peak-valley difference is 40.2%. Average natural inflow and total storage energy of hydropower plants are 264 m 3 /s and 4.65 × 10 6 MWh, respectively.

Effects of Inferring the Initial Values of Lagrange Multipliers
Conventionally, initial values of multipliers are set to the marginal cost corresponding to the solution of a simplified economic dispatch, or set to zero directly. Therefore, to test the effect of initial multipliers, three techniques are compared: (1) Zero value (ZV), where multipliers are set to zero directly; (2) A simplified version of economic dispatch (SED), which relaxes the integer constraints and is described in [15]; (3) The proposed method by expert system (ES1).
Moreover, to keep the single variable principle, PBM is adopted to update multipliers during the process. Table 2 shows the comparison of different initial multiplier generation techniques. In the stage of initial multiplier generation, the initial dual value of ES1 is 103,000 t, far more than those of ZV and SED, which means that ES1 provides the best lower bound of dual value at the beginning period. In the stage of iterative calculations, dual values of three techniques converge to a close value, around 111,303 t. Figure 5 shows the evolution process of dual value by three techniques. In Figure 5, as the iteration proceeds, the difference in dual value caused by initial multipliers becomes smaller and smaller. However, the number of iterations differs a lot, causing the differences in computational time, which are 782 s, 423 s, and 359 s for ZV, SED, and ES1, respectively. Moreover, the primal values of the three techniques are basically the same. The primal value of ES1 is only 0.03% less than that of ZV.  Therefore, the conclusions are drawn: (1) the quality of initial values of multipliers has a significant influence on the computational time, but little influence on the final objective function value; and (2) initial values of Lagrange multipliers obtained by ES1 are better than those by ZV and SED, demonstrating the superiority of ES1. Therefore, the conclusions are drawn: (1) the quality of initial values of multipliers has a significant influence on the computational time, but little influence on the final objective function value; and (2) initial values of Lagrange multipliers obtained by ES1 are better than those by ZV and SED, demonstrating the superiority of ES1.

Comparison with Standard PBM in Different Generation Scenarios
To demonstrate the robustness of IPBM, 12 typical generation scenarios from 2019 are chosen, of which the load demand and water inflow characteristics are summarized in Table 5. The scenarios are distinguished as follows: the first number refers to the scenario index (from 1 to 12), the second term is associated with the load demand level (LLD = low load demand; MLD = medium load demand; HLD = high load demand), and the third term means the water inflow condition (DS = dry season; WS = wet season). The conclusion is obtained by comparing the performance of IPBM with PBM in 12 scenarios.  In summary, compared with SG, CP, and PBM, ES2 shows extraordinary stability during the iteration process. Due to the advantage of the expert system, ES2 obtains the best objective function value in the least time, demonstrating its superiority in global search ability and computational efficiency.

Comparison with Standard PBM in Different Generation Scenarios
To demonstrate the robustness of IPBM, 12 typical generation scenarios from 2019 are chosen, of which the load demand and water inflow characteristics are summarized in Table 5. The scenarios are distinguished as follows: the first number refers to the scenario index (from 1 to 12), the second term is associated with the load demand level (LLD = low load demand; MLD = medium load demand; HLD = high load demand), and the third term means the water inflow condition (DS = dry season; WS = wet season). The conclusion is obtained by comparing the performance of IPBM with PBM in 12 scenarios. Table 6 lists the results for STHS by PBM and IPBM. From Table 6, it can be noticed that IPBM outperforms PBM in all scenarios in terms of dual value, primal value, and computational time. For dual value, the maximum and minimum relative differences are 0.64% in scenario 5_LLD_DS and 0.01% in scenario 4_LLD_DS, respectively, indicating that IPBM can provide a better lower bound for a primal problem. In terms of primal value, the average relative difference in 12 scenarios is −0.89%, which means that the average coal consumption reduction is 748 t by IPBM. Moreover, the computational time of IPBM is less than that of PBM in each scenario, verifying its high computational efficiency.   Figure 8 illustrates the relative difference between the primal value and natural inflow in different scenarios. It can be observed that the relative difference is approximately positively correlated with the natural inflow. From scenario 2_LLD_DS to scenario 5_LLD_DS, both the natural inflow and relative difference are at a low level. From scenario 6_MLD_WS to scenario 10_LLD_WS, the natural inflow and relative difference are at a high level, and gradually decrease. The maximum relative difference is 1.36% in scenario 6_MLD_WS. The relative difference in the primal value reflects the improvement in the objective function of IPBM compared with PBM. Therefore, the improvement in the IPBM objective function is more significant in large-inflow scenarios. The reason is that when the natural inflow is high, the corresponding hydropower potential is large. IPBM can improve the water use efficiency and yield more hydropower energy. As a result, coal consumptions in large-inflow scenarios are greatly reduced. at a high level, and gradually decrease. The maximum relative difference is 1.36% in scenario 6_MLD_WS. The relative difference in the primal value reflects the improvement in the objective function of IPBM compared with PBM. Therefore, the improvement in the IPBM objective function is more significant in large-inflow scenarios. The reason is that when the natural inflow is high, the corresponding hydropower potential is large. IPBM can improve the water use efficiency and yield more hydropower energy. As a result, coal consumptions in large-inflow scenarios are greatly reduced.

Conclusions
Hydrothermal power accounts for a high proportion in power systems, making STHS a challenge for operators. To efficiently solve the STHS problem, an IPBM combining the ES technique and standard PBM within the LR framework is presented in this paper. In IPBM, the ES consists of two parts: a knowledge base and inference engine. The knowledge base is built by extracting knowledge expressions from historical generation scenarios. Based on the knowledge base, initial values of Lagrange multipliers and new updates of Lagrange multipliers when null steps occur are reasoned by the inference engine. In this way, the ES provides a good lower bound to the primal problem at the beginning period. Moreover, the ES avoids solving quadratic programs when a null step is declaimed and increases the likelihood of reaching global optimum. To verify the effectiveness of IPBM, a case study is conducted in a large-scale hydrothermal system in China. Results in different cases indicate that, compared with several common methods, IPBM can obtain generation schedules with less coal consumption in less time, demonstrating its superiority in global search ability, computational efficiency, and robustness. Hence, IPBM proves to be an effective method for STHS.
Due to the negative impact of dam construction on river ecosystems, the ecological dispatch of hydropower is attracting increasing attention from researchers. Therefore, for further studies, it is necessary to bring ecological flow constraints into the present solution framework.  Data Availability Statement: Some data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request, including the code of the proposed IPBM, MILP, k-means, and silhouette coefficient.

Conflicts of Interest:
The authors declare no conflict of interest. Table 1 describes the basic parameters of thermal units. Table 2 lists the basic parameters of hydropower plants.