Solving the Short-Term Scheduling Problem of Hydrothermal Systems via Lagrangian Relaxation and Augmented Lagrangian

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 compatiblemodelling 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.


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 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 on/off 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.

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 min The objective function 2.1 seeks to minimize the total operating cost, which is composed by the immediate cost from stage 1 to T plus the expected future cost from the stage T 1 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 T 1 on.This cost is modeled by a multivariate piecewise linear function, which is dependent on the final volumes of the reservoirs at the stage T , as detailed ahead.

Hydraulic Constraints (C H )
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 , phjrt • 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 ph jrt q jrt , Q rt , s rt .

2.7
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 C H above in the compact form C H C HH α, Q, s, V ∩ C HUC z, q, Q, s, P H .The vectors z, q, Q, s, PH, and V gather the respective variables.The set C HH • represents constraints given by 2.2 -2.4 , while C HUC • represents 2.5 -2.8 .The objective function is given by f T pt, u, x f H α , where u and pt are vectors gathering all binary and continuous thermal variables, respectively.

Thermal Constraints C T
Let us consider the following equations: Constraints 2.9 describe the output limits for each unit.

2.10
Equation 2.10 represent the minimum uptime and downtime constraints

2.11
Equation 2.11 represent the ramp constraints.In a compact formulation, constraints in the set C T correspond to C T pt, u, x .

Hydrothermal Constraints C HT
We have the following equations:

2.12
These constraints represent the balance between supply and demand, for each time stage, and subsystem

2.13
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 C HT is written as C HT pt, PH,Int , where the vector Int gathers the exchanges variables.

Solution Strategy
Phase 1 Decomposition by Lagrangian Relaxation .The STS problem 2.1 -2.13 is given by We introduce artificial variables pta and PHa, which duplicate pt and PH, respectively.pta and PHa are used in constraints C HT • to substitute pt and PH.Variables Qa and sa duplicate Q and s, respectively.Qa and sa replace Q and s in pt pta, P H PHa, Q Qa, s sa.

3.2
Now, the artificial constraints that keep coupling are relaxed by introducing Lagrange

3.4
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 r and stage t 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 notglobally 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 s.t.: C T pt, u, x ∩ C HT pta, P Ha, Int ∩ C HH α, Qa, sa, V ∩ C HUC z, q, Q, s, P H .

3.10
Above c 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.

3.16
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 y is obtained in the following way: 3.17 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 x k is a vector with the solution of all variables for each k.

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.227GW.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 c 0 is equal to 1 × 10 −4 and it is increased along the iterations ξ as follows: where β 1 1.5, β 2 200, and β 3 10 4 .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 < β 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 c ξ 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 g x ξ 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.

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 .

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 PHD PHa − PH and QD Qa − Q at the last iteration in specific hydroplant with 3,300 MW installed capacity.The PHD bigger values are 2.39 MW, in stage 2. On the other hand, all QD values are null.

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 PHa − PH and Qa − Q at the last iteration in the same hydroplant

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%.

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.

Mathematical Problems in Engineering
The dual function A.3 can be rewritten as The maximum value of θ LR is 72.67 and λ 1 λ 2 36.33.The primal variables associated are pt 1 pt 2 3, pta 1 2, pta 2 0, and u 1 u 2 1. Clearly we have an infeasible primal solution.By using the AL shown in this paper, we build up the following AL dual function: A.5 The next step is to employ the APP technique to break the coupling of the AL dual function: A.6 The constants k 1 and k 2 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 k 1 and k 2 are 2.5 and 1.5, respectively.The AL dual function A.6 can be evaluated by means of three separable subproblems: where A.8 We use a gradient method to maximize the AL dual function: We set the initial value for ρ k as 1/50, and we update this penalty parameter by means of ρ k 1 2 ρ k .Proceeding in this way, the decomposition strategy presented in our paper converges in 37 iterations and finds the following primal-dual solution: pt 1 pta 1 2, pt 2 pta 2 u 2 0, u 1 1, λ 1 −15.29, and λ 2 36.40.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.

Constants and Indexes
i: Index of thermal units, so that r 1, I I: Number of thermal plants r: Index of reservoirs, so that r 1, R R:

Number of hydro plants t:
Index of stages, so that t 1, T T : Number of stages h a 0i , a 1i , a 2i : Fuel cost coefficients of the thermal plant i R$ Brazilian Real, as it is known the currency in Brazil., R$/MW and R$/MW 2 b 0i , b 1i : Cold fixed startup cost of the thermal plant i R$ ω i : Thermal time constant of the thermal plant i h A: Conversion factor of m 3 /s in units hm 3 y rt : Incremental inflow in the reservoir r and stage t m 3 /s v r min max : Minimum max.volume of the reservoir r hm 3 s r max : Maximum spillage of the reservoir r m 3 /s m: Index of reservoir upstream of the reservoir r τ mr : Number of stages that the outflow in the upstream hydro m takes to reach the downstream hydro r r : Set of reservoirs immediately upstream to the reservoir r P : Number of linear constraints of the future cost function FCF p: Index associated with the FCF π r p : Angular coefficient of the FCF R$/hm 3 α 0 p : Linear coefficient of the FCF R$ J r : Total of hydrounits associated with the reservoir r j: Index of hydrounits Φjr: Number of nonforbidden zones of the unit j and reservoir r k: Index The number of stages that the thermal unit has been on or off until the stage t h v rt : Volume of the reservoir r in the beginning of the stage t hm 3 q jrt : Turbined outflow of the unit j, reservoir r in the stage t m 3 /s Q rt : Turbined outflow in the reservoir r in the stage t m 3 /s s rt : Spillage of the reservoir r in the stage t m 3 /s ph jrt • : Power of the hydrounit j, reservoir r in the stage t MW Int elt − Int elt D et .

Figure 7 :
Figure 7: Different demand scenarios and primal recovery performance.
of nonforbidden zones of the hydrounits ph jkrt min max : Minimum maximum power of the hydrounit j, zone k, reservoir r and stage t MW pt i min max : Minimum maximum power of the thermal plant i MW rh rt : Power reserve of the hydroplant r and stage t MW t i up down : Minimum time required to start up shut down the thermal plant i h δ i : Ramp down rate of the thermal plant i MW Δ i : Ramp up rate of the thermal plant i MW e: Index of subsystems I e , R e : Set of all thermal/hydroplants of the subsystem e Ω e : Set of subsystems interconnected with subsystem e Int let max : Maximum power exchange from subsystem e to l in the stage t MW D et : Load of the subsystem and stage t MW .Binary variable which indicates if thermal unit i is operating u it 1 or not u it 0 during the stage t x it :