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

Authors: Rafael N. Rodrigues; Edson L. da Silva; Erlon C. Finardi; Fabricio Y. K. Takigawa
Journal: Mathematical Problems in Engineering
 (2012)

---

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


---

## Body

## 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 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. Sections2 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(2.1)min⁡FSTS=∑t=1T∑i=1Iuit[(a0i+a1iptit+a2iptit2)+(boi(1-exi,t-1/ωi)+b1i)(1-ui,t-1)]+α.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.


### 2.1. Hydraulic Constraints (CH)

Consider the following equations:(2.2)vr,t+1-vrt+A(Qrt+srt-∑m∈R+(r)(Qm,t-τmr+sm,t-τmr)-yrt)=0,(2.3)vrmin⁡≤vrt≤vrmax⁡,0≤srt≤srmax⁡,(2.4)α+∑r=1Rπr(p)vr,T+1≥α0(p),(2.5)Qrt-∑j=1J(r)qjrt=0.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,(2.6)∑k=1Φjrphjkrtmin⁡zjkrt≤phjrt(qjrt,Qrt,srt)≤∑k=1Φjrphjkrtmax⁡zjkrt.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](2.7)(∑j=1J(r)phj1rtmax⁡-PHrt)≥rhrt,PHrt=∑j=1J(r)phjrt(qjrt,Qrt,srt).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(2.8)zjkrt∈{0,1},∑k=1Φjrzjkrt≤1.Finally, (2.8) represents the integrality constraints.We write constraintsCH above in the compact form CH=CHH(α,Q,s,V)∩CHUC(z,q,Q,s,PH). The vectors z, q, Q, s, PH, and V gather the respective variables. The set CHH(·) represents constraints given by (2.2)–(2.4), while CHUC(·) represents (2.5)–(2.8). The objective function is given by fT(pt,u,x) + fH(α), where u and pt are vectors gathering all binary and continuous thermal variables, respectively.


### 2.2. Thermal Constraints(CT)

Let us consider the following equations:(2.9)ptimin⁡uit≤ptit≤ptimax⁡uitConstraints (2.9) describe the output limits for each unit.(2.10)uit={1if1≤xit<tiup,0if-1≥xit>-tidown0or1otherwise,,xit={max⁡(xi,t-1,0)+1ifuit=1,min⁡(xi,t-1,0)-1ifuit=0.Equation (2.10) represent the minimum uptime and downtime constraints(2.11)δi(ui,t-1,xit)≤ptit-pti,t-1≤Δi(ui,t-1,xit).Equation (2.11) represent the ramp constraints. In a compact formulation, constraints in the set (CT) correspond toCT(pt,u,x).


### 2.3. Hydrothermal Constraints(CHT)

We have the following equations:(2.12)∑i∈Ieptit+∑r∈RePHrt+∑l∈Ωe(Intelt-Intelt)=Det.These constraints represent the balance between supply and demand, for each time stage, and subsystem(2.13)0≤In⁡tlet≤In⁡tletmax⁡,0≤In⁡telt≤Inteltmax⁡.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 (CHT) is written as CHT = (pt, PH,Int), where the vector Int gathers the exchanges variables.


## 2.1. Hydraulic Constraints (CH)

Consider the following equations:(2.2)vr,t+1-vrt+A(Qrt+srt-∑m∈R+(r)(Qm,t-τmr+sm,t-τmr)-yrt)=0,(2.3)vrmin⁡≤vrt≤vrmax⁡,0≤srt≤srmax⁡,(2.4)α+∑r=1Rπr(p)vr,T+1≥α0(p),(2.5)Qrt-∑j=1J(r)qjrt=0.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,(2.6)∑k=1Φjrphjkrtmin⁡zjkrt≤phjrt(qjrt,Qrt,srt)≤∑k=1Φjrphjkrtmax⁡zjkrt.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](2.7)(∑j=1J(r)phj1rtmax⁡-PHrt)≥rhrt,PHrt=∑j=1J(r)phjrt(qjrt,Qrt,srt).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(2.8)zjkrt∈{0,1},∑k=1Φjrzjkrt≤1.Finally, (2.8) represents the integrality constraints.We write constraintsCH above in the compact form CH=CHH(α,Q,s,V)∩CHUC(z,q,Q,s,PH). The vectors z, q, Q, s, PH, and V gather the respective variables. The set CHH(·) represents constraints given by (2.2)–(2.4), while CHUC(·) represents (2.5)–(2.8). The objective function is given by fT(pt,u,x) + fH(α), where u and pt are vectors gathering all binary and continuous thermal variables, respectively.


## 2.2. Thermal Constraints(CT)

Let us consider the following equations:(2.9)ptimin⁡uit≤ptit≤ptimax⁡uitConstraints (2.9) describe the output limits for each unit.(2.10)uit={1if1≤xit<tiup,0if-1≥xit>-tidown0or1otherwise,,xit={max⁡(xi,t-1,0)+1ifuit=1,min⁡(xi,t-1,0)-1ifuit=0.Equation (2.10) represent the minimum uptime and downtime constraints(2.11)δi(ui,t-1,xit)≤ptit-pti,t-1≤Δi(ui,t-1,xit).Equation (2.11) represent the ramp constraints. In a compact formulation, constraints in the set (CT) correspond toCT(pt,u,x).


## 2.3. Hydrothermal Constraints(CHT)

We have the following equations:(2.12)∑i∈Ieptit+∑r∈RePHrt+∑l∈Ωe(Intelt-Intelt)=Det.These constraints represent the balance between supply and demand, for each time stage, and subsystem(2.13)0≤In⁡tlet≤In⁡tletmax⁡,0≤In⁡telt≤Inteltmax⁡.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 (CHT) is written as CHT = (pt, PH,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
(3.1)min⁡fT(pt,u,x)+fH(α),s.t.:CT(pt,u,x)∩CHT(pt,PH,Int)∩CHH(α,Q,s,V)∩CHUC(z,q,Q,s,PH).
We introduce artificial variables pta and PHa, which duplicate pt and PH, respectively. pta and PHa are used in constraints CHT(·) to substitute pt and PH. Variables Qa and sa duplicate Q and s, respectively. Qa and sa replace Q and s in CHH(·)(3.2)min⁡fT(pt,u,x)+fH(α),s.t.:CT(pt,u,x)∩CHT(pta,PHa,Int)∩CHH(α,Qa,sa,V)∩CHUC(z,q,Q,s,PH),pt=pta,PH=PHa,Q=Qa,s=sa.
Now, the artificial constraints that keep coupling are relaxed by introducing Lagrange multipliers λPT, λPH, λQ, λS:
(3.3)max⁡Θ(λPT,λPH,λQ,λS),(3.4)Θ(⋅)=min⁡fT(pt,u,x)+fH(α)+〈λPT,pt-pta〉+〈λHT,PH-PHa〉+〈λQ,Q-Qa〉+〈λs,s-sa〉,s.t.:CT(pt,u,x)∩CHT(pta,PHa,Int)∩CHH(Qa,sa,V)∩CHUC(z,q,Q,s,PH).
Above 〈·,·〉 denotes the dot product. The dual function (3.4) can be split as the sum of four separable terms, as follows:
(3.5)Θ(⋅):=θT+θHT+θHH+θHUC,(3.6)θT=min⁡fT(pt,u,x)+〈λPT,pt〉,s.t.:CT(pt,u,x),(3.7)θHT=min⁡-〈λPT,pta〉-〈λPH,PHa〉,s.t.:CHT(pta,PHa,Int),(3.8)θHH=min⁡fH(α)-〈λQ,Qa〉-〈λs,sa〉,s.t.:CHH(α,Qa,sa,V),(3.9)θHUC=min⁡〈λPH,PH〉+〈λQ,Q〉+〈λs,s〉,s.t.:CHUC(z,q,Q,s,PH).
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 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
(3.10)Φ:=min⁡fT(pt,u,x)+fH(α)+〈λPT,pt-pta〉+〈λHT,PH-PHa〉+〈λQ,Q-Qa〉+〈λs,s-sa〉+c‖pt-pta‖2+c‖PH-PHa‖2+c‖Q-Qa‖2+c‖s-sa‖2,s.t.:CT(pt,u,x)∩CHT(pta,PHa,Int)∩CHH(α,Qa,sa,V)∩CHUC(z,q,Q,s,PH).
Abovec 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:
(3.11)‖pt-pta‖2≈‖pt-ptaξ+ptξ2‖2+‖ptaξ+ptξ2-pta‖2,
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:
(3.12)max⁡Φ:=ϕT+ϕHT+ϕHH+ϕHUC,(3.13)ϕT=min⁡fT(pt,u,x)+〈λPT,pt〉+c‖pt-ptaξ+ptξ2‖2,s.t.:CT(pt,u,x),(3.14)ϕHT=min⁡-〈λPT,pta〉-〈λPH,PHa〉+c‖ptaξ+ptξ2-pta‖2+c‖PHaξ+PHξ2-PHa‖2,s.t.:CHT(pta,PHa,Int),(3.15)ϕHH=min⁡fH(α)-〈λQ,Qa〉-〈λs,sa〉+c‖Qaξ+Qξ2-Qa‖2+c‖saξ+sξ2-sa‖2,s.t.:CHH(α,Qa,sa,V),(3.16)ϕHUC=min⁡〈λPH,PH〉+〈λQ,Q〉+〈λs,s〉+c‖Q-Qaξ+Qξ2‖2+c‖s-saξ+sξ2‖2+c‖PH-PHaξ+PHξ2‖2,s.t.:CHUC(z,q,Q,s,PH).
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)y=∑kδkxk.
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 xk 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.Table 1
Thermal data.
PlantFuelFuel costsStart-up costsPower limitsMin. timesRampa0a1a2b0b1ωptminptmaxtuptdownδΔ01Nuclear5,624.020.90.000245,00012,5004706579915018002Nuclear5,624.016.30.000256,00021,00061251,3509920025003Gas4,000.0153.80.0005420250120092322703504Gas1,200.0117.00.00092851501806701110025005Gas657.074.40.005285421908691210030006Gas723.011.80.00439861604801112020007Coal986.0170.50.00893,2261,61371213246507008Coal1,183.0129.00.00773,8951,94582526258607509Coal1,540.075.00.00624,6912,3459201001020505010Gas755.3143.00.006588501656382215025011Gas529.070.20.005161351716382215025012Gas866.080.60.0066974812334711100200In the PR phase, the initial value ofc0 is equal to 1×10-4 and it is increased along the iterations ξ as follows:(4.1)cξ+1={β1cξifcξ≥β3c0,cξ+β2c0ifcξ≥β3c0,
where β1 = 1.5, β2 = 200, and β3 = 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 < β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(4.2)λξ+1=λξ+φ⋅c⋅g(xξ),
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.


### 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.Figure 1
LR phase convergence.We show the difference between the demand and the total generation for each stage in Figure2. 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).Figure 2
Infeasibility in demand attainment.


### 4.2. Primal Recovery

Figure3 shows the AL performance in the PR phase.Figure 3
PR phase convergence.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%.Figure4 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 4
Feasibility in demand attainment.Figure5 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.Figure 5
Hydrogeneration and turbined outflow infeasibility-PR phase.


### 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 (PHa-PH) and (Qa-Q) 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.Figure 6
Hydrogeneration and turbined Outflows infeasibility.


### 4.4. Sensitivity Analysis for Other Operation Conditions

To accomplish this task, we now consider four different demand scenarios. In Figure7, +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.Figure 7
Different demand scenarios and primal recovery performance.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 (Figure4). In all cases presented above, the value of the duality gap estimative, on average, was approximately 4%.


## 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.Figure 1
LR phase convergence.We show the difference between the demand and the total generation for each stage in Figure2. 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).Figure 2
Infeasibility in demand attainment.


## 4.2. Primal Recovery

Figure3 shows the AL performance in the PR phase.Figure 3
PR phase convergence.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%.Figure4 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 4
Feasibility in demand attainment.Figure5 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.Figure 5
Hydrogeneration and turbined outflow infeasibility-PR phase.


## 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 (PHa-PH) and (Qa-Q) 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.Figure 6
Hydrogeneration and turbined Outflows infeasibility.


## 4.4. Sensitivity Analysis for Other Operation Conditions

To accomplish this task, we now consider four different demand scenarios. In Figure7, +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.Figure 7
Different demand scenarios and primal recovery performance.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 (Figure4). 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.



---

*Source: 856178-2012-02-07.xml*