Abstract
This paper addresses the short-term scheduling problem of hydrothermal power systems, which results in a large-scale mixed-integer nonlinear programming problem. The objective consists in minimizing the operation cost over a two-day horizon with a one-hour time resolution. To solve this difficult problem, a Lagrangian Relaxation (LR) based on variable splitting is designed where the resulting dual problem is solved by a Bundle method. Given that the LR usually fails to find a feasible solution, we use an inexact Augmented Lagrangian method to improve the quality of the solution supplied by the LR. We assess our approach by using a real-life hydrothermal configuration extracted from the Brazilian power system, proving the conceptual and practical feasibility of the proposed algorithm. In summary, the main contributions of this paper are (i) a detailed and compatible modelling for this problem is presented; (ii) in order to solve efficiently the entire problem, a suitable decomposition strategy is presented. As a result of these contributions, the proposed model is able to find practical solutions with moderate computational burden, which is absolutely necessary in the modern power industry.
1. Introduction
The operation coordination of a hydrothermal system requires the use of several models given that a single model cannot accommodate all the complexities involved in this problem. Usually the problem is divided into scheduling problems with different modeling and planning horizons. The long-term scheduling model results are used as input for medium-term model which, in turn, feed their results to the short-term scheduling model. In the Short-Term Scheduling (STS) problem, a detailed modeling is necessary so that practical solutions can be found. To address this problem, a large-scale optimization problem with a mix of discrete and continuous variables needs to be formulated.
Some primal methods have been proposed for the STS based on, for example, progressive optimality [1], network flow [2], Benders decomposition [3], and mixed integer programming [4–6]. However, considering the detailed model used in our problem formulation and considering its constraints structure, dual decomposition appears as an approach for the problem. In this sense, Lagrangian Relaxation (LR) [7–19] and Augmented Lagrangian (AL) [20–22] have been extensively used. The main advantages of the LR are as follows: the original problem can be split into a sequence of smaller easy-to-solve subproblems, and a lower-bound for the optimal objective function is supplied. However, there is an important disadvantage: for nonconvex problems, the LR fails to find a feasible solution. On the other hand, with AL, it is possible to obtain a feasible solution or a near-feasible solution. However, the problem is not separable and it is necessary to use some inexact variant of this methodology.
Although almost all works which use dual decomposition in the STS problem are based on either LR [7–18] or AL [20–22], in this work we propose a two-phase approach similar to [23, 24] but contrasting to these references we introduced contributions, which are presented at end of this section. In the first phase, we use the LR to obtain (infeasible) primal solution, and in the second phase we use the AL to obtain a solution whose quality can be assessed by the LR. In [23], a two-phase approach is used, where the problem is decomposed into thermal unit commitment and linear hydrothermal scheduling subproblems. This decomposition is not efficient when the hydroproduction function is modeled by a nonlinear function and hydro-mixed-integer constraints are included in the STS problem. Aiming to deal with nonlinearities and binary variables of the hydrosystem, in our previous works [18, 19], a suitable decomposition was presented. In this paper, we use the LR of [18, 19] but including the AL second-phase optimization in order to achieve primal feasibility, modeling all nonlinearities and binary variables related to hydro- and thermal units.
Typically, LR technique relaxes coupling constraints such as demand and reserve requirements. Nevertheless, using this strategy, the hydro-subproblem remains coupled in time and space. An alternative approach consists in combining LR with Variable Splitting-LRVS method [25, 26], where the decomposition is achieved by duplicating some variables. In this paper, the LRVS is used to duplicate thermal and hydrovariables, as well as the turbined outflow and spillage variables. Thus, we obtain four subproblems: thermal, hydro, hydrothermal, and hydraulic. The first two subproblems take into account the unit commitment constraints (thermal and hydro). The hydrothermal subproblem considers demand reserve and transmission constraints. Finally, all the reservoir constraints are modelled into the hydraulic subproblem.
Given that the LRVS fails to find a feasible solution, an AL approach is used in attempt to overcome this issue. The artificial constraints relaxed in the LRVS are taken into account in the AL function. In order to maintain the decomposition, the Auxiliary Problem Principle (APP) [27] is employed. We assess our approach by using a real-life hydrothermal configuration extracted from the Brazilian power system, proving the conceptual and practical feasibility of the proposed algorithm.
The main contributions of this paper are the following.
(i) Although [18, 19] present a detailed modelling for this problem, these works did not solve the entire STS problem. Actually, in the works [18, 19], only the hydroelectric subproblem was solved and considering only the LR phase. It was shown explicitly the necessity of a robust primal recovery phase since the primal solution supplied by the LR cannot be applied in real-life decisions. Furthermore, the presented results shown in [18, 19] did not take into account the Future Cost Function, which makes the problem more difficult to be solved since all reservoirs, even in different cascades, are mathematically coupled. In this paper, we solved the entire problem, which includes the thermal unit commitment subproblem and constraints to satisfy the demand.
(ii) Regarding [23], which proposes a two-phase approach to solve the entire problem, it is worth to point out that in this paper the problem related to the hydroelectric system is modelled as continuous linear problem, which can lead to nonpractical solutions, given that the power production function is clearly nonlinear. Furthermore, it is necessary to include 0-1 variables in the problem to take into account the onoff unit’s status as well as the forbidden zones. This issue is very important for a hydro-dominated system and all these issues are considered in this paper.
(iii) Furthermore, in terms of solution strategy, differently of [18, 19, 23, 24], in this paper we include the pseudoprimal point strategy [28] as a warm starting by the APP. This strategy improves dramatically the performance of the APP algorithm.
Given that the problem under consideration is classified as a large-scale one and it is very complex due to several nonlinearities associated with hydro- and thermal technologies, this work focuses on dual decomposition algorithms. Alternative techniques, based on Mixed Integer Linear Programming (MILP), have been used in other applications, especially for mid-sized problems with thermal predominance or small-size hydroproblems [6, 29–32].
This paper is organized as follows. Sections 2 and 3 present the mathematical formulation and the solution strategy, respectively. Section 4 reports the numerical results and Section 5 provides some concluding remarks. The nomenclature is represented in the end of paper.
2. Problem Modeling
The problem to be solved in real life occurs in the context of a very tense operational process. Usually, in the power industry, the latest data collection phase ends close to 10:00 AM and a feasible schedule have to be computed by the Independent System Operator (SO) as soon as possible. On the other hand the solution of this problem needs to be unambiguous, which requires a detailed modelling so that postoptimization adjustments, involving human expertise, should be avoided. Thus, based on these requirements we present in the following a detail model for our problem; therefore, as consequence a suitable strategy solution is necessary aiming to reduce the computational burden.
The objective function for the STS problem is
The objective function (2.1) seeks to minimize the total operating cost, which is composed by the immediate cost from stage 1 to plus the expected future cost from the stage on. The immediate cost is composed by the variable, start-up, and fixed thermal costs [14]. On the other hand, the last term in (2.1), supplied by a medium-term model [33], is the expected future cost from the stage on. This cost is modeled by a multivariate piecewise linear function, which is dependent on the final volumes of the reservoirs at the stage , as detailed ahead.
2.1. Hydraulic Constraints ()
Consider the following equations:
Constraints (2.2), (2.3), (2.4), and (2.5) represent the stream-flow balance equations, the volume and spillage limits, the Future Cost Function (FCF), and the penstock water balance, respectively,
Constraints (2.6) describe the output limits for the hydro-units nonforbidden zones. In (2.6), is a high-order nonconvex polynomial, which allows to model hydraulic losses, nonlinear tailrace levels, and nonlinear turbine-generator efficiencies. More details about this function can be found in [18, 19]
Constraints (2.7) model the hydroreserve requirement and the reservoir power balance. Typically, in hydrothermal systems dominated by hydroplants, these plants provide spinning reserve better than thermal units. Hydrounits respond faster than thermal units to system disturbs given that they can increase or decrease their outputs quickly. Nevertheless, if the physical nature of the power system requires also thermal spinning, we do not see difficulty to include reserve constraints which consider additionally the thermal plants
Finally, (2.8) represents the integrality constraints.
We write constraints above in the compact form ∩. The vectors , , , , , and gather the respective variables. The set represents constraints given by (2.2)–(2.4), while represents (2.5)–(2.8). The objective function is given by (, , ) + (α), where and are vectors gathering all binary and continuous thermal variables, respectively.
2.2. Thermal Constraints
Let us consider the following equations:
Constraints (2.9) describe the output limits for each unit.
Equation (2.10) represent the minimum uptime and downtime constraints
Equation (2.11) represent the ramp constraints. In a compact formulation, constraints in the set correspond to.
2.3. Hydrothermal Constraints
We have the following equations:
These constraints represent the balance between supply and demand, for each time stage, and subsystem
Constraints (2.13) model the subsystems power exchange limits. This type of simplification is adopted in real life given the huge size of the Brazilian transmission system and its corresponding complexity, and they are supplied by power flow and transitory stability studies. The set () is written as = (, ,Int), where the vector Int gathers the exchanges variables.
3. Solution Strategy
Phase 1 (Decomposition by Lagrangian Relaxation). The STS problem (2.1)–(2.13) is given by We introduce artificial variables and , which duplicate and , respectively. and are used in constraints to substitute and . Variables and duplicate and , respectively. and replace and in Now, the artificial constraints that keep coupling are relaxed by introducing Lagrange multipliers , , , : Above denotes the dot product. The dual function (3.4) can be split as the sum of four separable terms, as follows: Subproblem (3.6) is a nonlinear mixed 0-1 optimization problem, which corresponds to the thermal unit commitment, and is solved by a dynamic programming algorithm similar to described in [34]. In this algorithm, the generation level was discretized in steps of 1% and this strategy has shown coherency with the practical expectation. However, this approach can lead to suboptimal results as shown in [35], and for different systems it is interesting to use other solution strategies to solve the thermal unit commitment subproblems such as presented in [35].Subproblems (3.7) and (3.8) are Linear Programming (LP) problems, which can be solved by any LP solver. Subproblem (3.9) is a nonlinear mixed 0-1 optimization problem, and it corresponds to the hydrounit commitment. In (3.9), all feasible unit combinations of a hydro and stage are solved aiming to find an optimal solution. To solve the continuous problems, we use a sequential quadratic programming (SQP) algorithm. Additional details of this algorithm can be seen in [18, 19]. Since our contributions are not focused on the maximization of the dual function (3.4), we use the N1CV2 Bundle subroutine [36] whose main algorithm does not take into account the disaggregated version of Bundle method which, as it can be seen in [37], can improve the LR performance. Additionally, the SQP algorithm used to solve (3.9) may supply not-globally solutions and, as a consequence, the dual function and subgradient vector values can be evaluated inexactly. In this sense, there are more efficient Bundle methods which can deal with this inaccuracy, as it can be seen in [38] and references therein.
Phase 2 (Primal Recovery by Augmented). In order to search a primal feasible solution, we use an inexact Augmented Lagrangian-AL, similar to [23, 39]. Starting from (3.2), the same artificial constraints are relaxed, and the new dual function is presented below Above is the penalty parameter. Notice that the quadratic term makes the dual function nonseparable. By applying the APP [27], the following approximations are carried out: where the variables with ξ correspond to the values obtained at the previous iteration. Similarly, the same strategy (3.11) must be done for the other quadratic terms in (3.10). Then, the dual problem related to AL can be rewritten as follows: It can be observed that the subproblems (3.13)–(3.16) present the same structure as in the LR phase, but with an important difference; now the LPs (3.7) and (3.8) are transformed into quadratic programming problems (3.14) and (3.15) (In the Appendix, we show how in detail the proposed two-phase approach works.). Given that the APP approach is valid only locally, it is essential to have a good primal starting point available. In this paper, we use the pseudoprimal point found by convex combination of the Bundle actives cuts [28]. In general terms, the pseudoprimal point is obtained in the following way: In (3.17), k is the number of Bundle actives cut, δ is a constant, between 0 and 1 of the convex combination (the sum of all values of vector δ is 1), and is a vector with the solution of all variables for each k.
4. Numerical Results
We assess the solution strategy on a real-life hydrothermal configuration extracted from the Brazilian power system. More precisely, we consider a system with 121 hydro- and 12 thermal units whose maximum installed capacity is 38.227 GW. Five cascades make up the system; the largest one possesses seven plants with 58 units and the smallest one is composed of two plants with 12 units. Detailed data for the hydrosystem are too lengthy to list in this paper and can be found in [19, 40]. We show the data for the thermal plants in Table 1. The two-day planning horizon is discretized hourly. Initial reservoir volumes were taken at 50% of usable volumes. The interconnected hydrothermal system is divided into four subsystems.
In the PR phase, the initial value of is equal to and it is increased along the iterations ξ as follows: where = 1.5, = 200, and = 104. An important practical question is how to select the penalty parameter sequence. In this sense, the main considerations for selecting the parameter sequence are the following [41]: (i) 1 < ≤ 2 when in the AL maximization there are hard-to-solve primal subproblems; (ii) the penalty parameter should be increased exponentially in the initial iterations and, in the sequence, it is interesting to have a less ambitious update scheme. Given that the primal solution supplied by the LR possesses a large number of infeasible constraints, the use of an exponential rule of in the beginning is recommended. This approach speeds up the convergence process. The algorithm may switch to a linear update when the gradient vector norm reaches some tolerance.
When used the steps (i) and (ii), even changing the parameters, important variations in the CPU time were not observed. On the other, for different data conditions (for instance, for different demand and inflows scenarios), a fast and easy adjustment of the constants in (3.17) may be necessary.
In this paper, the multipliers are updated as where φ is the step size and is the constraint value.
This strategy avoids the inherent oscillatory behavior of the gradient method over the iterations. The algorithm stopping criteria are defined by an absolute error tolerance ε between each original variable and its associated duplicated variable.
4.1. Optimization by LR
The initial Lagrange multipliers are λPT = λPH = λQ = λS = 1.0. We found a primal function value of R$ 1,124,688,100.00 after 450 iterations, which took 92 minutes of CPU time in an Intel Core 2 Extreme CPU X9650 @ 3.00 GHz, 4.00 GB of RAM memory. In Figure 1, the subgradient vector norm and the dual function over the iterations are depicted. Although the dual function is close to the optimal value by iteration 300, the subgradient norm is still dropping in further iterations.
We show the difference between the demand and the total generation for each stage in Figure 2. Note that the deviations are not so small. For instance, in stage 25 the deviation between demand and total generation is approximately 2.5 GW. These differences cannot be simply adjusted in real-time operation. A great part of this infeasibility is related to the oscillatory aspect in the LR when applied to LP such as subproblems (3.7) and (3.8).
4.2. Primal Recovery
Figure 3 shows the AL performance in the PR phase.
The algorithm stopping criteria are defined by ε = 2%, that is, each gradient vector component must be less than 2% of the maximum value of the associated variable. After 103 iterations and 11 minutes of CPU time, the operating cost obtained was R$ 1,193,547,056.00. Since it is not possible to evaluate precisely the duality gap (we do not have the global optimal solution), we can compute the gap between the value of the dual function computed in Phase 1 and the total cost of the primal solution obtained in Phase 2. From weak duality, this gap gives a bound for the optimal cost. In this first case, the relative gap is 13.89%.
Figure 4 shows the difference between the demand and the total generation for each hour. As it can be seen, the algorithm supplies a highprecision feasible solution.
Figure 5 shows the values of and at the last iteration in specific hydroplant with 3,300 MW installed capacity. The bigger values are 2.39 MW, in stage 2. On the other hand, all values are null.
4.3. Sensitivity Analysis
A new stopping criterion is used to carry out a sensitivity analysis of the AL performance, where ε = 10%. We found an optimal value after 47 iterations that took 5 minutes. Figure 6 shows the values of and at the last iteration in the same hydroplant considered in Figure 5. The bigger values are 19.73 MW, in stage 1, and 0.005 m3/s, in stage 20, respectively.
4.4. Sensitivity Analysis for Other Operation Conditions
To accomplish this task, we now consider four different demand scenarios. In Figure 7, +10%, −10%, and −20% represent, the demand of case base, presented in Figure 4, plus 10%, minus 10%, and minus 20%, respectively. In the tests, the same rules (3.17), (4.1), and (4.2) were used.
As it can be noticed, regardless of the operation, point the algorithm shows be convergent in all scenarios, without any change in their rules and parameters. Specifically in these cases, the computational burden did not increase significantly, in comparison with the base case (Figure 4). In all cases presented above, the value of the duality gap estimative, on average, was approximately 4%.
5. Conclusions
The STS problem formulation involving simultaneous LR and AL techniques are presented in this paper. LR is used to decompose the problem into four subproblems whose scheme is able to address properly the nonlinear production function and mixed 0-1 constraints of the hydrounits. Since the primal problem is nonconvex, we apply an inexact AL aiming to find a feasible primal solution. The infeasibility between generation and demand is eliminated and the artificial constraints are satisfied. The authors are still investigating other forms to update the penalty parameter. The objective is to improve the solutions and the measure of their quality. However, the results in this paper are promising for the application of this methodology for huge problems such as the entire Brazilian power system.
Appendix
The problem presented in [42] is used to show how our solution strategy works. We rewrote the original problem as
Using the splitting variable technique, the problem is written as follows: Relaxing the artificial constraints in (A.2): The dual function (A.3) can be rewritten as The maximum value of θLR is 72.67 and = = 36.33. The primal variables associated are = = 3, = 2, = 0, and = = 1. Clearly we have an infeasible primal solution. By using the AL shown in this paper, we build up the following AL dual function: The next step is to employ the (APP) technique to break the coupling of the AL dual function: The constants and are calculated as the variables average value obtained in the previous iteration. In the first iteration, the AL uses the primal variables obtained in the last iteration of the LR. In consequence, the initial values for and are 2.5 and 1.5, respectively. The AL dual function (A.6) can be evaluated by means of three separable subproblems: where We use a gradient method to maximize the AL dual function: We set the initial value for as 1/50, and we update this penalty parameter by means of . Proceeding in this way, the decomposition strategy presented in our paper converges in 37 iterations and finds the following primal-dual solution: , , , and . The corresponding objective function value is 104. Indeed, this is the optimal solution of the original problem. To summarize the convergence process, Figure 8 shows the dual function and the gradient vector norm over the iterations.
List of Symbols
Constants and Indexes| : | Index of thermal units, so that , |
| : | Number of thermal plants |
| : | Index of reservoirs, so that , |
| : | Number of hydro plants |
| : | Index of stages, so that , |
| : | Number of stages (h) |
| , , : | Fuel cost coefficients of the thermal plant (R$ (Brazilian Real, as it is known the currency in Brazil.), R$/MW and R$/MW2) |
| , (): | Cold (fixed) startup cost of the thermal plant (R$) |
| : | Thermal time constant of the thermal plant (h) |
| : | Conversion factor of (m3/s) in units (hm3) |
| : | Incremental inflow in the reservoir and stage (m3/s) |
| : | Minimum (max.) volume of the reservoir [hm3] |
| : | Maximum spillage of the reservoir (m3/s) |
| : | Index of reservoir upstream of the reservoir |
| : | Number of stages that the outflow in the upstream hydro takes to reach the downstream hydro r |
| : | Set of reservoirs immediately upstream to the reservoir |
| : | Number of linear constraints of the future cost function (FCF) |
| : | Index associated with the FCF |
| : | Angular coefficient of the FCF (R$/hm3) |
| : | Linear coefficient of the FCF (R$) |
| : | Total of hydrounits associated with the reservoir |
| : | Index of hydrounits |
| Φ: | Number of nonforbidden zones of the unit and reservoir |
| : | Index of nonforbidden zones of the hydrounits |
| : | Minimum (maximum) power of the hydrounit , zone , reservoir and stage (MW) |
| : | Minimum (maximum) power of the thermal plant (MW) |
| : | Power reserve of the hydroplant and stage (MW) |
| : | Minimum time required to start up (shut down) the thermal plant (h) |
| : | Ramp down rate of the thermal plant (MW) |
| : | Ramp up rate of the thermal plant (MW) |
| : | Index of subsystems |
| , : | Set of all thermal/hydroplants of the subsystem |
| : | Set of subsystems interconnected with subsystem |
| : | Maximum power exchange from subsystem to in the stage (MW) |
| : | Load of the subsystem and stage (MW). |
| α: | Expected future cost value ($) |
| : | Binary variable which indicates if thermal unit is operating () or not () during the stage |
| : | The number of stages that the thermal unit has been on or off until the stage (h) |
| : | Volume of the reservoir in the beginning of the stage (hm3) |
| : | Turbined outflow of the unit , reservoir in the stage (m3/s) |
| : | Turbined outflow in the reservoir in the stage (m3/s) |
| : | Spillage of the reservoir in the stage (m3/s) |
| : | Power of the hydrounit , reservoir in the stage (MW) |
| : | Power of the hydroplant in the stage (MW) |
| : | Binary control variable which indicates if the hydrounit of the reservoir is operating () or not () in the zone during the stage |
| : | Power of the thermal plant in the stage (MW) |
| : | Power exchange from subsystem to in the stage (MW). |