Solving the unit commitment problem of hydropower plants via Lagrangian Relaxation and Sequential Quadratic Programming

 

 

Erlon C. Finardi^I; Edson L. da Silva^I; Claudia Sagastizábal^II

^ILabPlan - Laboratório de Planejamento de Sistemas de Energia Elétrica Campus Universitario, Trindade, Florianópolis - 88040-900, SC, Brasil, E-mail: erlon@labplan.ufsc.br/ edson@labplan.ufsc.br
^IIIMPA - Instituto de Matematica Pura e Aplicada Estrada Dona Castorina 110, Jardim Botânico, Rio de Janeiro - 22460-320, RJ, Brasil, E-mail: sagastiz@impa.br

 

 

---

ABSTRACT

We consider the optimal scheduling of hydropower plants in a hydrothermal interconnected system. This problem, of outmost importance for large-scale power systems with a high proportion of hydraulic generation, requires a detailed description of the so-called hydro unit production function. In our model, we relate the amount of generated hydropower to nonlinear tailrace levels; we also take into account hydraulic losses, turbine-generator efficiencies, as well as multiple 0-1 states associated with forbidden operation zones. Forbidden zones are crucial to avoid nasty phenomena such as mechanical vibrations in the turbine, cavitation, and low efficiency levels. The minimization of operating costs subject to such detailed constraints results in a large-scale mixed-integer nonlinear programming problem. By means of Lagrangian Relaxation, the original problem is split into a sequence of smaller and easy-to-solve subproblems, coordinated by a dual master program. In order to deal better with the combinatorial aspect introduced by the forbidden zones, we derive three different decomposition strategies, applicable to various configurations of hydro plants (with few or many units, which can be identical or different). We use a Sequential Quadratic Programming algorithm to solve nonlinear subproblems. We assess our approach on a real-life hydroelectric configuration extracted from the south sub region of the Brazilian hydrothermal power system.
Mathematical subject classification: 90C90, 46N10, 47N10.

Key words: hydrothermal systems, Unit Commitment Problems, Lagrangian Relaxation, sequential quadratic programming.

---

 

 

1 Introduction

The optimal generation scheduling is an important daily activity for electric power generation companies. The goal is to determine which units are to be used in order to generate enough power to satisfy demand requirements and various technological constraints, with minimum operating cost. In particular, hydrothermal systems must consider the stream-flow equations for reservoirs. These equations couple all the reservoir along a hydro-valley, because the amount of outflow water¹ released by one power plant affects water volumes in all the plants downstream. Furthermore, water travel times and alternative uses of water, such as irrigation or flood control, for example, must also be taken into account.

The optimal scheduling of hydropower plants is called the Hydro Unit Commitment (HUC) problem. To solve the HUC problem, a highly sophisticated modeling for the operation of hydro plants is required. Specifically, a hydropower plant may be composed of several turbine-generator groups, referred to in this work as ''units''. The amount of power generated by one hydro unit depends on the efficiency of both the turbine and the generator, as well as on the net head and the unit turbined outflow. In turn, the net head is a nonlinear function of the storage and of the reservoir outflow. The joint turbine-generator efficiency varies with the net water head and the unit turbined outflow. In addition, the existence of forbidden operation regions prevents the unit from generating power in a wide and continuous range. These regions, modeled by 0-1 variables, aim at avoiding vibrating modes that may produce unwanted power oscillations, cavitation phenomena, and low levels of efficiency. Thermal power plants have simpler production functions, but they need start-up and shut-down times and they often present nonlinear operating costs.

As a result, the hydrothermal unit commitment problem is a large-scale mixed-integer nonlinear programming problem which can only be effectively solved by applying decomposition techniques. Lagrangian Relaxation (LR) is particularly suitable for this type of problems, [1-5], although some other methodologies have also been proposed [6]. However, so far none of the works in the area has considered a modeling as comprehensive as ours, with representation of hydraulic losses, nonlinear tailrace levels, turbine-generator efficiencies, and forbidden operation zones.

With respect to the solution method, our contribution consists in a thorough analysis of three different decomposition schemes, all derived from LR. The first strategy relies on a complete enumeration of all possible 0-1 operating states of the units composing a hydropower plant. This approach is suitable for plants with a low number of identical units. The second strategy, requiring less computation effort, is applicable for plants with many units that are different and have many forbidden regions. The third strategy combines the two other approaches, and can be used in systems with both types of hydro plants.

Many of the subproblems resulting from our decomposition schemes are Nonlinear Programs (NLP) of small size. We solve them by a Sequential Quadratic Programming (SQP) method [7,8], in a quasi-Newton variant [9,10], which presents good convergence properties.

Our work is organized as follows. In Section 2 we give the hydro plants and units modeling. Section 3 is devoted to the mathematical formulation of the optimal scheduling problem. The solution strategy, with the three decomposition schemes, is given in Section 4. In Section 5 we report numerical results on a hydrothermal system corresponding to Brazil's southern electric sub region. We end in Section 6 with some concluding remarks.

 

2 Hydro generating units

For a unit j, the generated power, ph_j, expressed in [MW], depends on the unit turbined outflow, q_j, and on the net water head, hl_j, and on joint turbine and generator efficiency, h_j:

The net water head has the expression:

Here, fcm stands for the forebay level. For short term horizon problems, as in our case, the forebay level remains practically constant, specially in the Brazilian case, whose huge reservoirs have typical regularization levels of a couple of years. Therefore, we consider fcm constant. By contrast, the downstream level fcj(.) varies abruptly in short times, mainly due to the plant turbined outflow, Q, given by the addition of the outflows of all the units composing the plant. For some power plant configurations, fcj also varies with the reservoir spillage, s. In (2), the term k_j represents hydraulic losses resulting from friction of the water in penstock, where k_j is a constant expressed in s²/m⁵ [11].

The unit efficiency, depending on hl_j and q_j, is usually represented by hill diagrams given by the factory; see Figure 1. We estimate it by interpolation; see [12], using a polynomial function:

where the coefficients r_0j,..., r_5j have been computed beforehand.

 

 

Figure 1 also displays some important operating constraints on the turbine-generator group. For example, for net head values smaller than the so-called nominal level (41,5m), the turbine is unable to make the generator attain its nominal power (120 MW). On the other hand, for values higher than 41,5m, there is a limit of power limit imposed by the generator capabilities, because the turbine could effectively reach power levels beyond 120 MW. Since at some power levels cavitation phenomena and nasty mechanic vibrations may appear, in order to extend the lifetime of the unit and to avoid power oscillations, such power levels are forbidden. For example, Figure 1 shows a forbidden operation region ranging from 70 to 90 MW.

By combining (1)-(3) we obtain our model for the hydro production function:

For the Brazilian case, fcj(Q,s) is represented by a fourth degree polynomial. Therefore, from (4) we see that ph_j is a polynomial² of degree 12 on the variables Q and s, and of order 7 in the variable q_j.

At first sight our model may appear as ''too complicated''; however, it is important to realize that only such a detailed description can accurately represent the diverse amounts of power generated by a unit at different operating states.

 

3 Problem formulation

The objective function for the thermal-HUC problem has the expression:

Here, the planning horizon is composed by T time steps, the thermal mix has I plants, c_it(.) represents the operating cost of the i-th thermal plant at time step t, and a stands for the system expected future cost at the end of the planning horizon; see (8) below. Frequently c_it(.) includes fixed costs as well as fuel costs related to start-up and nominal generation of thermal units [2],[5].

We formulate the thermal-HUC constraints by splitting them into three different subsets, C_H, C_T and C_HT, corresponding to the respective variables involved namely, hydraulic, thermal, or both. Each subset is characterized by a specific type of coupling, such as units in the same power plant along different time steps (time coupling), or different power plants in a given time step (space coupling).

We now proceed to give each constraint in detail.

3.1 Constraints involving only hydraulic variables (C_H)

• Stream-flow balance equation:

We use the index r for reservoirs, v is the reservoir storage, y is the incremental inflow, is a set gathering all reservoirs upstream the r-th, and t_mr is the water travel time between reservoirs m and r.

• Maximum and minimum storage, and maximum spillage per reservoir:

• Expected future cost function, given by longer term planning models, and estimating the cost of using today water that might become necessary (and expensive) in the future; see [13]. It is a piecewise affine function that depends on the final levels of stocked water, v_rT:

• Penstock water balance equation per reservoir:

J(r) is the number of generating units in reservoir r.

• Power limits, given for each operating region of the unit:

F _jr denotes the total number of non-forbidden regions of the j-th unit in reservoir r; k is the corresponding index, and stand for the minimum and maximum power limits. The binary variable z_jkrt is 1 if the j-th unit in reservoir r is operating in the k-th region at time step t, and it is set to 0 otherwise.

• Reservoir power balance:

• Reserve constraints:

rh_rt is the minimum reserve of reservoir r at time step t.

• Integrality constraints:

In the sequel, to alleviate notation, we write constraints C_H above in the abstract form

where the vectors z, q, Q, s, PH and V gather the respective variables. The set C_HH represents constraints given by (6)-(8), modeling the reservoirs, while C_HUC represents the unit constraints, i.e., (9)-(13). In this abstract formulation, a = a(V).

3.2 Constraints involving only thermal variables (C_T)

• Power limit for each unit:

Here stand for the minimum and maximum power limits of unit i. The binary variable u_it is 1 if the unit is operating at time step t, and it is set to 0 otherwise.

• Reserve constraints:

rt_rt is the reserve of unit i at time step t.

• Minimum up-time, ti^up, and downtime, ti^down, for each unit:

where the state variable x_it is equal to the number of time steps the unit has been up/down until time t.

• Ramp constraints:

d _i(u_i,t-1 ,x_it) < pt_it -pt_i,t-1 < D_i (u_i,t-1 ,x_it) ,

• d _i(.) and D_i(.) are the maximum allowed variations of generation of the unit between two time steps.

In an abstract formulation, constraints in the set C_T correspond to C_T(u,pt), where u and pt are vectors gathering all binary and continuous thermal variables, respectively.

3.3 Constraints involving both hydraulic and thermal variables (C_HT)

• Satisfaction of demand, per time step and subsystem:

The interconnected hydrothermal system is divided into subsystems, indexed by e. Accordingly, all thermal units (reservoirs) of subsystem e are gathered in the index set I_e (R_e). There are W_e subsystems interconnected with subsystem e; the exchange of energy at time t, is denoted by Int_let (Int_elt respectively) when it goes from subsystem l to e (from e to l, respectively). Finally, D_et is the demand of subsystem and at time t.

• Subsystems exchange limits, from e(l) to l(e), at time t, , ():

In our abstract notation, the set C_HT is written as C_HT (pt,PH,Int), where the vector Int gathers the subsystem exchanges.

The above description confirms the level of complexity of the optimization problem to be solved. We now address the solution strategy adopted in this work.

 

4 Solving the HUC problem

The economic impact of the optimal scheduling of power plants is undeniable. Because of their solid theoretical background, LR techniques appear in this area as the preferred solution method. In particular, multipliers associated to demand constraints given by (14) are used to price energy.

The divide to conquer approach of LR, also called price decomposition [9], is well known. Essentially, coupling constraints are relaxed via Lagrange multipliers whose corresponding dual problem is decomposable into simpler subproblems (called local subproblems). The coordination of subproblems is then done by a master program, which finds new multipliers by making one iteration of a nonsmoth algorithm that maximizes the dual function.

There are many ways of relaxing coupling constraints. An important criterion for deciding how to proceed is the resulting duality gap, which should be the smallest possible. In this matter, the introduction of artificial variables to uncouple constraints appears as a good choice; see [14], and also [5,15] for an application to the thermal UC problem. For this reason, we apply a similar approach in this work, and derive three different decomposition schemes, adapted to different unit configurations in the Brazilian hydrothermal system.

4.1 First decomposition strategy - D1

In the abstract notation, the thermal HUC problem becomes:

To achieve decomposition, we introduce artificial variables pta and PHa, which duplicate, respectively, pt and PH. Variables pta and PHa are used in constraints C_HT to replace pt and PH. In addition, artificial variables Qa and sa duplicate Q and s, respectively. Qa and sa replace Q and s in C_HH. With these additional variables, (15) is rewritten as follows:

In (16) the newly introduced artificial constraints hold the coupling of the problem. Hence, we relax them by associating Lagrange multipliers l_PT, l_PH, l_Q, l_S and writing the corresponding dual problem³:

Problem (17) can be rewritten as follows:

where:

In the LR approach, the primal problem (15) is replaced by the dual problem (18) whose objective function D1(l), can be split as the sum of four terms, corresponding to subproblems (19)-(22). Subproblem (19) is a nonlinear optimization problem with continuous and binary variables, coupled along time steps, but not along plants. It can be solved by a classic Dynamic Programming method; as in [1],[5]. Subproblem (20) is a standard linear programming (LP) problem, coupled along plants, but not along times steps, which can be solved by any (LP) commercial solver. Subproblem (21) is also an LP problem, coupled both in time and space via the stream-flow constraints given by (6). Even though (21) can be large-scale, an LP solver can still solve it efficiently. Finally, subproblem (22) is a nonlinear mixed-integer optimization problem, uncoupled both in time and units. This subproblem corresponds to the commitment of hydro units, for a given reservoir and time step. The higher J(r) and F_jr (the number of units in the reservoir and of operating zones, respectively) are, the bigger computational effort will be required to solve (22).

Each sub-subproblem in (22), for each time step and for a given power plant, is a mixed-integer NLP problem, with binary variables corresponding to different operating modes in the plant. The total number of possible operating modes is given by the product of all combinations of the operating modes of all the units composing the plant. Each combination of a unit is a configuration where the corresponding binary variables are fixed to one of the feasible values. Once the binary values are fixed, the problem becomes a nonlinear program, whose size is dependent on J(r).

Generally, hydropower plants have identical units, and each unit has a single operating zone. In this case, the total number of modes is no longer 2^J(r), but J(r)+1 and, thus, a complete enumeration of modes seems a good strategy. Sometimes, however, there are power plants with many different types of units, and several operating modes. For these configurations, an enumeration procedure may become too expensive from the computational point of view. We now introduce an alternative decomposition scheme, adapted to such situations.

4.2 Second decomposition strategy - D2

In order to avoid the enumerative process required to solve subproblem (22), we eliminate the coupling between the binary variable z and the continuous variables [q,Q,s] which appear in C_HUC. Therefore, we rewrite (15) as follows:

Now the set C_HUCa gathers constraints given by (9) and (11), C_HUCb contains (10) and (13) and C_HUCres corresponds to the reserve constraint (12). Besides the artificial variables used in (16), we use pha to replace the hydro production function ph(q,Q,s) in the set C_HUCb, and rewrite (23) as:  

Now not only the artificial constraints, but also the constraint set C_HUCres keep the problem coupled. Hence, we relax these constraints by introducing multipliers l_PT,l_PH, l_Q, l_S, l_ph and l_Res and writing the dual problem:

Here, ph^max and rh are vectors corresponding to a unit's maximum capacity and to reserve levels for the plant, respectively. Separability in (24) has now the expression:

where D1_T(l), D1_HT(l), D1_HH(l) are the dual functions from (19), (20) and (21). The remaining terms are:

Subproblem (25) is an NLP problem (with only continuous variables), for each time step and power plant. Subproblem (26) is a mixed-integer linear program on variables corresponding to a single unit, which can be solved by enumeration of the operating zones.

4.3 Third decomposition strategy - D3

Our last decomposition combines D1 and D2. Accordingly, we employ D1 for those plants with a reduced number of units and operating zones, while D2 is applied for plants for which D1 would be too expensive. We split the hydro plants index set into R(G) and R(X), corresponding, respectively, to plants where (22) and (25)-(26) are applied:  

We proceed like for D1 and D2, but splitting the reservoir index sets:

After relaxation, the dual problem of (27) is:

The dual functions D1_T(l), D1_HT(l), D1_HH(l) are those in (19), (20) and (21), respectively. The dual function D1_HUC(l) from (22) only applies for reservoirs with index r Î R(G). The remaining dual functions in (28) apply to reservoirs r Î R(X), and are given by subproblems (25)-(26).

4.4 Nonlinear programming subproblems

A crucial issue for an effective application of LR is the fast resolution of subproblems defining the dual function for a fixed multiplier l. Since many of such subproblems are nonlinear programs, we implemented a Sequential Quadratic Programming (SQP) quasi-Newton method. More precisely, each iteration k of the algorithm generates a direction p^k by solving the quadratic programming problem:

Here, f(x^k), c_e(x^k) and c_i(x^k) represent, respectively, the objective and equality and inequality constraint functions at a point x^k. The matrix M^k estimates the Lagrangian Hessian for the NLP, L^k. In order to avoid the calculation of second-order derivatives, and to preserve positive definiteness of the sequence of quasi-Newton matrices, we use a BFGS (Broyden-Fletcher-Goldfarb-Shanno) [16], formula, appended with a Powell correction [17]:

where:

Globalization of the method is achieved by performing a line search on the following function⁴:

known as Han's merit function [9]. Since f(x^k,s^k) is an exact penalty function, there is a finite positive value such that an unconstrained minimum of f(x^k,s^k) solves the original NLP for all s^k > . We update the parameter s^k accordingly. More precisely, the directional derivative of the merit function along the direction satisfies the relation:

where z is the Lagrange multiplier associated with constraints in (29). For any stationary point of (30), such as p^k, it can be shown that the estimate in (31) gives a descent direction for f if M^k is positive definite and s^k is updated in order to satisfy:

We exit the line search when the Armijo [18] condition is satisfied:

here, w Î ]0,1/2[ and a is the positive stepsize. Ideally, D^k should be the exact value of the directional derivative. We estimate it by the upper bound from (31):

D ^k : = Ñf(x^k)^T p^k- s^k||c( x^k)^#||_¥

Finally, we add an extra term to p^k in order to avoid Maratos effect [19]. This phenomenon may impair the superlinear local convergence rate by rejecting unit stepsizes when close to a solution. The corrected direction p^k, as shown in [9], ensures that asymptotically there is enough constraint reduction; see 5.1 below for more details.

 

5 Numerical Results

We assess the three decomposition schemes on a real-life hydroelectric configuration extracted from the Southern region of the Brazilian hydrothermal power system. More precisely, we consider a system with 121 generating units whose maximum installed capacity is 31.129,2 MW. ⁵ Figure 2 reports the data for the system, where power plants numbered #3, #4, #6 and #14 have production functions independent of spillage. Values between brackets in Fig. 2 correspond to water travel times, expressed in hours. For this configuration, the biggest power plants are #7 and #16, with 20 units each one, while the smallest plant, #14, has only 2 units. The planning horizon of two days is discretized in hourly time steps, yielding T = 48. Initial reservoir volumes were taken at 50% of the usable volumes, while the inflows were considered null. We do not address here the dual solution in detail, we refer to [5] for this subject. Instead, we fixed Lagrange multipliers for each reservoir and time step, l_PHrt, and use a proximal quasi-Newton variant of a bundle method to optimize the remaining multipliers; see Ch. 9 in [9]. The values for l_PHrt are chosen based on generation costs associated with typical demand curves, i.e., with higher values for peak times with high demand⁶.

 

 

We implemented the three dual subproblems D1_UCH, D2_UCH, D3_UCH, corresponding respectively to (22), (25)-(26) and (28). For the dual solution we use N1CV2 code; see [20]. Subproblem D1_UCH has (R+Rv)×T variables, where Rv denotes the number of reservoirs with production function depending on spillage. Since T = 48, R = 18 and Rv = 14, subproblem (22) has 1536 variables. In (25)-(26) subproblem D2_UCH has (2×R+Rv+ngmix)×T = 8208 variables, where ngmix = 121 denotes the total number of units in the mix. The size of subproblem D3_UCH depends on the decomposition scheme. The LP given by D1_HH (21), as well as other LPs are solved using ILOG CPLEX 7.1 solver. For D1_UCH there are, in addition, 6288 NLP problems. Finally, for D2_UCH there are R×T = 864 NLP and T×ngmix = 5808 easy mixed-integer LP problems ⁷.

For D1_UCH we found an optimal value of $ -17.339.230,0, after 175 iterations, that took 180 minutes of CPU times in a Pentium III 550 MHz computer with 128 Mb of RAM memory. For D2_UCH the optimal value found was $ -17.450.168,0, after 325 iterations in 50 minutes. Since the dual value in D2_UCH is smaller than the one from D1_UCH, and dual values give lower bounds for the primal optimal value, we can conclude that primal variables associated with D1_UCH are better than those associated with D2_UCH.

To further assess the previous remark, we now consider in more detail some selected primal variables, for some specific reservoirs. Figures 3 and 4 show the values for Q and Qa obtained at the last iteration of both D1_UCH and D2_UCH for Sobradinho #18 power plant. We also show the value of l_PHrt (price) used to solve the subproblems.

 

 

 

 

It can be seen in the figures that Qa takes mostly two values: 0 and 4278 m³/s (maximum value). This bang-bang behaviour is explained by the linear nature of subproblem (21). The values for Q exhibit a different behaviour, closer to the price profile, i.e., to l_PHrt. From the comparison of D1_UCH and D2_UCH we see that the relaxed primal constraint⁸ is more violated for D2_UCH, a problem that contains a higher number of relaxed constraints. However, infeasibility becomes smaller for time steps with bigger l_PHrt, i.e., for peaks of demand. For these time steps, primal points obtained with D2_UCH are good approximations to those from D1_UCH, requiring a computational effort that can be up to 3 times bigger. For lower prices, both problems give primal points that are even more infeasible, the worse values being associated with D2_UCH.

The results observed for Sobradinho power plant are typical for all the plants in the mix, with variations in the computed primal infeasibility. In general, we observed that for reservoirs downstream, that tend to operate with outflows near to the nominal values, primal points obtained from D2_UCH were close to those from D1_UCH. We report this behaviour in Figures 5 and 6, with the results for Ilha Solteira - 7 power plant.

 

 

 

 

We conclude from our experiments that, in terms of primal solutions, subproblem D1_UCH is a better option. However, since average hourly CPU timesfor D1_UCH are 3.5 higher than D2_UCH, the enumerative process required byD1_UCH should not be used for power plants with complex configurations, such as #7 and #16. In Table 1 we report the main results for the third decomposition scheme, D3_UCH, where for R X = {7,16} we applied the scheme corresponding to (25)-(26). We also give, for comparison purposes, the previously obtained values for D1_UCH and D2_UCH.

 

 

The computational effort required for solving subproblems involving the big plants #7 and #16 (each one having 20 identical units , is clear in Table 1. Even though D3_UCH¹⁰ needs to solve 27,50% NLP problems less than D1_UCH, in terms of CPU times the gain was of 61,11%. Another important matter shown by Table 1 concerns dual optimal values: the (absolute value) difference between D2_UCH and D1_UCH is of $ 110.938, while the difference between D3_UCH and D1_UCH is $ 19.847, i.e., about 5.5 smaller.

After solving the dual problem, it is necessary to adjust the primal points corresponding to the dual optimum in order to find a primal point that is feasible, i.e., that satisfies the constraints. In order to recover primal feasibility, a purification-like phase should be executed afterwards. Such processes are often based on heuristics depending on the particular problem structure; see for example [21,22]. In particular, for the Brazilian case, general-purpose combinatorial optimization heuristics are not suitable. An augmented Lagrangian technique seems in this case better. We refer to [5] for a description of this technique in a similar context.

5.1 Sequential Quadratic Programming Algorithm

In our work, NLP problems have ng+2 variables and 2ng+2 constraints, where ng is the number of generating units in the considered configuration (in our example ng < 20). SQP algorithms do not need starting feasible points; in our implementation starting points only satisfy constraints (6) and (10)(i.e., both the penstock stream-flow equations and power limit constraints), but not the reserve constraint (12). In addition, to prevent the starting quadratic program (29) to be infeasible, we introduced additional constraints on the direction obtained from (29), aimed at satisfying physical operating bounds. Such bounds are related to each unit maximum turbined outflow, as well as maximum spillage levels. We mention that more sophisticated alternatives would be possible, for instance considering active constraints as [23], or introducing slack variables in a trust-region SQP, as in [24].

To solve each quadratic program (29), we use the Fortran code PLCBAS, an active set QP solver described in [25]. Once the direction is computed, the algorithm checks validity of (32) and then performs, if needed, a correction strategy to avoid Maratos effect. A remedy proposed by Maratos in [19] is to move towards the feasible region by using a second order model of the constraints. A more simple strategy is to use the corrected step p^k as follows:

where is a right inverse for the Jacobian matrix whose rows are the transposed gradients of the constraint functions at x^k. The correction can be thought of as a constraint restoration step from x^k+p^k. It is called 'second-order correction' because c(x^k+p^k) = O(|p^k|²). Although this modification preserves fast convergence, it may no longer yield a descent direction for the merit function (30). For this reason, instead of searching the stepsize along the line x^k+ap^k, we use the arc

Since p^k is tangent to the arc at a = 0, Armijo condition (32), written with p^k replaced by the corrected step p^k, is eventually satisfied whenever p^k is a descent direction, and the ''arc'' search ends; see Prop. 15.7 in [9]. The additional number of constraint evaluation required by our correction is compensated by the robustness and efficiency gained by the method. In our runs, at each NLP solution, the second order correction was used in average 30% of the iterations.

We use the following stopping test:

opt_test = ||Ñ L^k||_¥ +||c(x^k)^#||_¥ < eps,

with eps = 1,0×10^-8. Emergency exits, after 300 iterations of 600 simulations were also implemented. The observed average values for achieving convergence were of 8 iterations and 10 simulations.

In order to assess our self-made algorithm, we compare its performances the NLP solver Easy!, available for academic applications from the Optimization Group at IME, University of Campinas, Brazil [26]. This solver uses the augmented Lagrangian method described in [27]. We obtained identical results, with inferior CPU times than those employed by Easy!. However, we did not take advantage of all resources available in Easy! solver¹¹. Furthermore, our own implementation was tailored for the structure of our specific problem.

 

6 Concluding Remarks

We address in this work the problem of optimal commitment of hydraulic generating units in a hydrothermal power system. Two main topics were discussed, namely the modeling and solution strategy for the hydraulic problem. We gave a detailed modeling for the production function of each hydro unit, which takes into account the effect of variable efficiency rates, hydraulic losses, tailrace levels, as well as multiple operating zones. With respect to the solution methodology, we gave an LR decomposition approach using variable duplication to uncouple difficult constraints. We assessed the decomposition method by implementing our approach and testing it on a real-life hydraulic configuration, extracted from the Brazilian system. Our implementation focused on the hydraulic subproblems. We analyzed in detail the practical applicability of three decomposition schemes, in terms of CPU times and obtained primal points. For our configuration, the enumerative process appeared preferable for most power plants. For the whole Brazilian system, double in size to our test, this approach might become too expensive and a combined strategy, like the one given in our third decomposition scheme, may be preferable.

Acknowledgments. Research supported by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and Centro de Pesquisas de Energia Elétrica (CEPEL). The authors acknowledge with grateful thanks the contributions made by Maria Elvira Maceira and Albert Cordeiro Geber de Melo from CEPEL.

 

REFERENCES

[1] S. Ruzic, N. Rajakovic and A. Vuckovic, A Flexible Approach to Short-Term Hydro-Thermal Coordination, PART I: Problem Formulation and General Solution Procedure, IEEE Transactions on Power Systems, 11(3), August (1996).         [ Links ]

[2] C. Li, E. Hsu, A.J. Svoboda, C.L. Tseng and R.B. Johnson, Hydro Unit Commitment in Hydro-Thermal Optimization, IEEE Transactions on Power Systems, 12(2), May (1997).         [ Links ]

[3] X. Guan, A. Svoboda and C.A Li, Scheduling Hydro Power Systems with Restricted Operating Zones and Discharges Ramping Constraints, IEEE Transactions on Power Systems, 14(1), February (1999).         [ Links ]

[4] E. Ni, X. Guan and R. Li, Scheduling Hydrothermal Power Systems with Cascaded and Head-Dependent Reservoirs, IEEE Transactions on Power Systems, 14(3), August (1999).         [ Links ]

[5] A. Belloni, A. Diniz, M.E.P. Maceira and C. Sagastizábal, Bundle Relaxation and Primal Recovery in Unit Commitment Problems - The Brazilian Case, Annals of Operations Research Special - Volume on OR models for Energy Policy, Planning and Management, Kluwer Academic Publishers, Netherlands, 2001.         [ Links ]

[6] G.B. Sheble and G.N. Fahd, Unit Commitment Literature Synopsis, IEEE Transactions on Power Systems, 9(1), February (1994).         [ Links ]

[7] P.T. Boggs and J.W. Tolle, Sequential Quadratic Programming, Acta Numerica, (1996), 1-51.         [ Links ]

[8] J. Nocedal and S.J. Wright, Numerical Optimization, Springer Series in Operations Research (1999).         [ Links ]

[9] J.F. Bonnans, J.C. Gilbert, C. Lemaréchal and C. Sagastizábal, Numerical Optimization. Theoretical and Practical Aspects, Universitext, Springer-Verlag, Berlin (2002).         [ Links ]

[10] D.P. Bertsekas, Nonlinear Programming, Athena Scientific, 2^nd Edition, Belmont, MA, (1999).         [ Links ]

[11] A. Arce, T. Ohishi and S. Soares, Optimal Dispatch of Generating Units of the Itaipú Hydroelectric Plant, IEEE Trans. Power Systems, 17(1), February (2002).         [ Links ]

[12] F.J. Heredia and N. Nabona, Optimum Short-Term Hydrothermal Scheduling with Spinning Reserve Through Network Flows, IEEE Trans. Power Systems, 10(3), August (1995).         [ Links ]

[13] E.L. da Silva and E.C. Finardi, Parallel Processing Applied To the Planning of Hydrothermal Systems. IEEE Transactions on Parallel and Distributed Systems, 14 (2003), 721-729.         [ Links ]

[14] C. Lemaréchal and A. Renaud, A Geometric Study of Duality Gaps, With Applications, Mathematical Programming, 90(3) (2001), 399-427.         [ Links ]

[15] C. Lemaréchal, F. Pellegrino, A. Renaud and C. Sagastizábal, Bundle Methods Applied to the Unit Commitment Problem, System Modeling and Optimization, (1996), 395-402.         [ Links ]

[16] S.P. Han, Superlinearly Convergent Variable Metric Algorithms for General Nonlinear Programming Problems, Mathematical Programming, 11 (1976), 263-282.         [ Links ]

[17] M.J.D. Powell, A Fast Algorithm for Nonlinearly Constrained Optimization Calculations, Numerical Analysis, 144-157, Springer-Verlag, Berlin (1978).         [ Links ]

[18] L. Armijo, Minimization of Functions Having Lipschitz Continuous First Partial Derivates, Pacific Journal of Mathematics, 16(1-3) (1966).         [ Links ]

[19] N. Maratos, Exact Penalty Function Algorithms for Finite Dimensional and Control Optimization, Ph.D. Thesis, University of London (1978).         [ Links ]

[20] C. Lemaréchal and C. Sagastizábal, Use of the Code N1CV2, Inria, France (2002).         [ Links ]

[21] H. Blondel and A. Renaud, Décomposition Par Les Prix et Lagrangien Augmenté: Une Association Efficace Pour le Traitement des Problèmes de Production, Bulletin EDF, Série C, 4 (1991), 23-43.         [ Links ]

[22] F. Zhuang and F.D. Galiana, Towards a More Rigorous and Practical Unit Commitment by Lagrangian Relaxation, IEEE Transactions on Power Systems, 3(2), May (1988).         [ Links ]

[23] C.T. Lawrence, A Computationally Efficient Feasible Sequential Quadratic Programming Algorithm, Ph.D. Thesis, University of Maryland, USA (1998).         [ Links ]

[24] M.R. Celis, J.E. Dennis and R.A. Tapia, A Trust-Region Strategy for Nonlinear Equality Constrained Optimization, Numerical Optimization, SIAM, Philadelphia, 71-82 (1985).         [ Links ]

[25] E. Casas and C. Pola, An algorithm for indefinite quadratic programming based on a partial Cholesky factorization, RAIRO-Recherche Opérationnelle/Operations Research, 27 (1993), 401-426.         [ Links ]

[26] Internet address: http://www.ime.unicamp.br/ ~ martinez/software.htm.

[27] A. Friedlander, J.M. Martinez and S.A. Santos, A New Trust Region Algorithm for Bound Constrained Minimization, Applied Mathematics and Optimization 30 (1994), 235-266.         [ Links ]

 

 

Received: 30/VIII/04. Accepted: 22/X/04.

 

 

#613/04.
Correspondence to: Erlon C. Finardi
1 Turbined and spilt water flow.
2 Since hl_j = f(,Q⁴,s⁴), it follows that h_j = f(,Q⁸,s⁸) and, thus, by (1), ph_j = f(,Q¹²,s¹²).
3 From now on, the Euclidean inner product of two vectors, l and v, will be denoted by l^Tv = S_i l_iv_i .
4 The symbol ^# is used to denote only active constraints at x^k.
5 This amount corresponds to about 49% of the total hydraulic capacity of Brazil.
6 These costs presented values between 0 and 45 $/MW.
7 These  are  indeed  easy  to  solve  problems,  with  only  two  variables  (one  integer,  one continuous),  and  two  constraints.
8 The difference between Q and Qa.
9 Total number of NLP problems solved at each iteration.
10 Ilha Solteira  and Tucuruí  6 power plants are the only ones where (22), requiring the enumerative process, was not employed; we use (25)-(26) instead.
11 Easy! Solver is more efficient when analytic derivatives are available for the simulations; otherwise a finite difference approximation needs to be estimated at each simulation.