Abstract
This paper presents a new approach for the expected cost-to-go functions modeling used in the stochastic dynamic programming (SDP) algorithm. The SDP technique is applied to the long-term operation planning of electrical power systems. Using state space discretization, the Convex Hull algorithm is used for constructing a series of hyperplanes that composes a convex set. These planes represent a piecewise linear approximation for the expected cost-to-go functions. The mean operational costs for using the proposed methodology were compared with those from the deterministic dual dynamic problem in a case study, considering a single inflow scenario. This sensitivity analysis shows the convergence of both methods and is used to determine the minimum discretization level. Additionally, the applicability of the proposed methodology for two hydroplants in a cascade is demonstrated. With proper adaptations, this work can be extended to a complete hydrothermal system.
1. Introduction
The long-term hydrothermal system operation (LTHSO) objective is to determine the amount each hydro- and thermal plant must generate for minimizing the expected operation cost at a given horizon of T monthly stages [1].
To optimize the use of water the hydroplants are constructed on the same river basin, allowing water inflows in one reservoir to be used in other plants for generating energy [2]. This causes a spatial dependency between reservoirs, as the operation of one hydroplant affects the availability of resources of the other reservoirs downstream [3, 4]. Since the decisions regarding water storage in one stage influences the availability of resources in future stages, this problem has already a time dependency constraint [5].
Moreover, the difficulty to forecast water inflows makes it impossible to plan exactly how much water should be stored in the reservoirs during rainy periods, as well as how much energy should be generated using thermal plants.
In large-scale systems the LTHSO problem becomes very complex due to nonlinearities, such as thermal costs or hydrogeneration functions, and stochasticity of state variables considering demands or water inflow [1]. A possible way to assess the impacts of present decisions on the future is to use an optimization algorithm, such as dynamic programming, to estimate the expected cost-to-go function.
To investigate this problem there were many studies in the 1970s and early 1980s, with the majority using dynamic programming [1, 4, 5]. For example, with the increasing Brazilian electrical system's complexity, the dynamic programming state variables grew exponentially, resulting in a phenomenon known as the “curse of dimensionality,” making it impossible to solve the real problem with reasonable accuracy. However, in 1985 the Stochastic Dual Dynamic Programming (SDDP), using Bender's decomposition [6], was proposed to deal with this curse of dimensionality [7, 8] solving a subsample of the complete problem as separate small ones. It is possible to model the state variables without entirely discretizing the state space, as the states used to model the cost-to-go function are evaluated by previous sample procedure.
Another technique for reducing the problem's dimension and complexity is based on the concept of aggregated reservoirs [9–11]. These reservoirs are an estimation of the energy generation of a group of hydroplants with similar water inflow characteristics. In summary, a group of hydroplants is modeled as a unique equivalent reservoir thereby drastically reducing the problem's dimension.
Other techniques to tackle the DP problem include the DP successive approximation (DPSA) method [12], where each reservoir is optimized at a time, assuming a fixed operation for the remaining reservoirs. The two-stage algorithm minimizes total costs from the first-stage decisions plus the total expected future cost, being the cost of all future decisions, which depends on the first-stage solution [13, 14]. Incremental Dynamic Programming and Differential Dynamic Programming [15] were also used in the reservoir optimization problem. These methods are generally useful techniques for the deterministic case; however they were not successful in the stochastic multireservoir case, as presented by Labadie [15].
Recently, there have been many advances in integrating the DP with other algorithms, mostly including heuristic techniques, such as neurodynamic programming [16, 17], genetic dynamic programming [18], and swarm optimization dynamic programming [19], with just a few applied to the LTHSO problem. For example, the GA was applied to the Brazilian hydrothermal system by Leite [20], producing significant results. Additionally, neural networks were used to model the cost-to-go functions with an efficient state space discretization scheme by Cervellera et al. [21], but using nonlinear optimization. The drawbacks of heuristic and nonlinear methodologies are that they require specialized system knowledge otherwise the optimization solution does not converge to an optimal solution [15].
Nowadays, the SDDP methodology is used in many countries, as in the case of the Brazilian power system, where the SDDP with aggregated reservoirs is still the official methodology used for determining the long-term hydrothermal system operation, the short-run marginal cost, among other applications. Nevertheless, alternatives must be sought to compare these results.
Although SDDP with the use of aggregated reservoirs can solve the problem in reasonable computer time, there is a price to be paid when the cost-to-go function is not well estimated for every important part of the problem's space state. The solution determined can be far from the optimal solution of a stochastic dynamic programming problem (SDP), where the space state is entirely discretized and used for estimating the expected cost-to-go function. In a single hydroplant study case, comparing SDP and SDDP, there were lower average hydroelectric generation and higher operation costs using the latter, as shown by Martinez and Soares [22].
There are two main ways used for solving dynamic programming problems for hydrothermal systems operation. The first uses the expected future costs table to coupling subproblems through stages. This approach is disadvantageous as the subproblems cannot be modeled as a linear programming problem (LP), and specialist algorithms have to be developed for each type of problem. The second one is to adjust hyperplanes for every bunch of points, allowing the subproblems to be modeled and solved as linear programming problems by using an LP solver. This approach has many advantages such as modeling simplicity, ability to solve large-scale problems, convergence to global optimal solutions, and using the capacity of modern LP solvers, which are developed with code optimization [15].
However, its main disadvantage is that it adds useless constraints to the problem by representing hyperplanes repeatedly. As the number of hyperplanes grows, the LP solver may take much time in getting to the optimal solution of the subproblems. Since it is necessary to solve many subproblems, any small inefficiency can result in a big time loss for completing execution of the dynamic programming problem.
In summary, this work proposes a new way to model the expected cost-to-go function for the DP using the Convex Hull (CH) algorithm [23, 24]. This approach makes it possible to use a table of expected future values as a piecewise linear function calculated by the CH algorithm, allowing the dynamic programming subproblems to be solved with LP solvers instead of some specialist algorithm. The cost-to-go function modeled by CH has the minimum number of hyperplanes needed to represent future cost-to-go table and, unlike the second approach described before, makes each LP subproblem easier to solve, resulting in less computational effort to solve the complete DP problem.
The outline of this paper is structured as follows. Section 2 presents a dynamic programming model applied to the long-term operation planning problem; Section 3 introduces the methodology, and Section 4 explains how the convex hull algorithm is used in the cost-to-go functions formation with a tutorial example. Section 5 illustrates some case studies. Conclusions are drawn in Section 6 and, finally, Section 7 describes possible future work.
2. Stochastic Dynamic Programming—Model Description
Dynamic Programming (DP) is a method for solving sequential decision problems, that is, complex problems that are split up into small problems, based on Bellman’s Principle of Optimality [25]. The hydrothermal operation planning problem is broken up into small subproblems, each one corresponding to each period. The optimal decision of the current stage depends on a series of future decisions. In other words, the future decisions depend strongly on the current decision with a direct impact on the total operation cost. Using the DP approach, the problem is solved by starting in the last decision stage with a recursion in time. The optimal solution in each stage balances the decision on that stage with future stages [26].
Long-term hydrothermal operational planning is, by nature, a stochastic problem as the incremental water inflows are uncertain. This leads to the Stochastic Dynamic Programming (SDP), a technique widely used in the LTHSO problem, as described by [15, 26–28]; see also [29].
Considering the state vector being the initial reservoir volume and the incremental water inflow in each hydroplant in the stage t, the SDP is formulated as the following optimization model [26, 30]: subject to where T is number of stages; expected operation cost of state incremental water inflow in stage (month); expected value considering the probability of water inflow in the stage conditioned to state state vector at the beginning of stage generation decision in the stage generation cost in the stage discount rate; thermoelectric generation (MW-month); hydroelectric generation (MW-month); load demand (MW-month); reservoir storage of hydroplant in the beginning of stage (hm3); reservoir storage at the end of stage (hm3); incremental water inflow in stage (hm3/month); water discharge in stage (hm3/month); maximum turbined outflow (hm3/month); water spillage in stage (hm3/month); maximum reservoir volume (hm3); minimum reservoir volume (hm3).
Expression (2.1) represents the objective function, where is the immediate operation cost of the decision This vector is composed of the reservoir volume at the end of the stage, turbined outflow, water spillage, as well as the thermal generation and unsupplied load costs. Additionally the objective function is composed of the future operation cost
The water balance is expressed in (2.3), where the volume stored at the end of a stage, is equal to the volume stored at the beginning of this stage, plus the water inflow to this reservoir during the stage, yt, less the discharge, and the amount of spilled water, This equation represents the coupling between successive stages.
The constraints (2.4) give the maximum and minimum reservoir storage levels, and the constraints (2.5) give the maximum outflow.
Consequently, the immediate operation costs are the thermal generation costs necessary to meet the load in stage Part of this demand uses hydrogeneration accordingly to the operation decision, and the remaining load is met using thermal generation.
This problem can be solved using DP, by the discretization of the state variables, analyzing all possible combination of those variables. This discretization is a linear piecewise approximation to the real cost-to-go function that presents nonlinearities due to thermal costs and hydroproduction functions. An example of this linear approximation is shown in Figure 1 for a single hydroplant system, where 5 and 21 discretizations are considered in order to model the cost-to-go function in a given month. The blue dots outline the optimal hydrothermal operation for that specific discretized reservoir level.
The disadvantage of this technique is the “curse of dimensionality,” being the exponential growth of the computational requirements with the increase in the dimension of the state vector. This problem reinforces the usage of new techniques to get the expected cost-to-go functions.
3. Proposed Methodology—Convex Hull
The Convex Hull (CH) algorithm calculates, given a finite set of points, the boundary of the minimal convex set containing those points [23]. For example, this algorithm was used by Diniz and Maceira [31] to improve modeling of the hydropower production function in the short-term hydrothermal dispatch problem of a large-scale system, resulting in a new four-dimensional piecewise linear model.
A set is said to be convex if, for all and every point connecting and is also in So is convex if and only if it contains all the convex combinations of its elements. That is,
A variety of algorithms are proposed for solving the convex hull, including [24] Graham scan, Jarvis march, Divide and Conquer, and QuickHull [23]. The latter is used in the development of this paper, where the Convex Hull methodology is used in order to model the expected cost-to-go functions, of the long-term hydrothermal system operation.
At each stage, the optimal dispatch is computed for each state space discretization, as a linear programming problem. The set of points for all the state space discretizations is used in the Convex Hull algorithm. This algorithm gives the set of hyperplanes forming a convex set, that will be used as a piecewise linear approximation for the expected cost-to-go functions at each stage of the stochastic dynamic programming problem. The algorithm steps can be seen on Figure 2.
Following the dynamic programming approach, the modeling of this function starts at the last stage Using linear programming, step 1 determines the mean optimal operation costs for each state. The set of points for the reservoir’s storage level with mean optimal cost is used in the convex hull algorithm for calculating the set of convex planes (or hyperplanes). These planes are used as a piecewise linear approximation for the expected cost-to-go functions (ECFs), as step 2 shows and, next, in step 3 there is a recursion in time. The planes in the previous step are used as optimization constraints in this iteration, and this procedure is repeated until it reaches the first stage, as shown in step 4. This proposed methodology is thoroughly explained with a single reservoir system example in the next section, Tutorial System.
4. Tutorial System
The computation of the expected cost-to-go function in the backward phase of the stochastic dynamic programming problem is used for analyzing the proposed methodology. A three-stage tutorial system, composed of 1 hydro- and 2 thermal power plants is used. Table 1 shows the hydropower plant characteristics and Table 2 the thermal plants characteristics [32].
The system load is 45 MW-month fixed for all stages, and the penalties for failure in load supply are represented by a dummy thermal plant, assumed as 1000 $/MW-month. Plus, two water inflow scenarios are considered as seen in Table 3.
In this example only 3 discretizations are considered for simplifying the analysis, which are 100%, 50%, and 0% of the useful reservoir capacity, as shown in Figure 3.
The first step consists of the computation of the mean optimal cost for each discretization level, being the mean of optimal costs calculated for each water inflow scenario, for each discretization, at each stage.
The stochastic dynamic programming approach starts with the last stage where there is no future cost. So the immediate cost, determined by the thermal generation and the penalty for the lost load, is the optimal cost calculated by an LP problem for each water inflow, as modeled in Section 2.
The amount of reservoir storage at the end of the stage, water discharge, thermal generation decision, and the amount of unsupplied load for each scenario, as well as the mean optimal cost for each discretization level, can be seen in Table 4. Notice the relationship between the amount of water in the reservoir and the water inflow with the immediate cost.
After calculating the mean optimal costs, these values are used in the Convex Hull algorithm to get a set of lines (planes or hyperplanes when a higher number of hydroplants are considered), which is a convex set, as seen in Figure 4.
By analyzing Figure 4, it is necessary to eliminate the line that links the extreme points with Figure 5 being a linear piecewise approximation to the expected cost-to-go function.
From this point, the approximated cost-to-go function of the stage 3 will be used as a future cost in the next step for calculating the expected cost-to-go function in stage 2. Due to nonlinearities of thermal costs and hydrogeneration functions, this is a linear piecewise approximation of the cost-to-go function [13, 22].
It is worth noting that as the convex hull algorithm delineates the boundary of the minimal convex set, repeated planes are not given. For example, in Figure 5, if the state space was also discretized for the 75% level of the reservoir, there is a null cost. Accordingly, two repeat lines would be given, one corresponding to the 50% and another one for the 75% level. The CH algorithm gives only one line for this set of 3 points, which represents the elimination of unnecessary redundant constraints in the LP problem.
The same is repeated for the second stage and is shown in Table 5, this time considering the future cost, with the respective expected-cost-to-go function seen in Figure 6.
With those expected cost-to-go functions, it is possible to calculate the expected cost for any discretization level in the first stage. This stage represents the current month, where the reservoir storage level is known. Keeping the same standard, Table 6 shows the mean optimal cost for the first stage, that represents the expected optimal cost for operating the hydrothermal system in the three stages, for the same discretization levels.
The approximate expected cost-to-go function for the whole horizon considered is shown in Figure 7. For example, if the reservoir storage level measured is 50%, the expected cost to operate this system in the three stages, considering the two scenarios, is $597.8.
5. Case Studies
In this section, a series of case studies composed of two cascaded reservoirs are proposed as shown in Figure 8.
5.1. System Configuration
Table 7 shows the hydropower plants characteristics used in the proposed study.
Additionally, 3 thermal plants are considered. Table 8 shows the main characteristics of these thermal plants. The minimum thermal generation is assumed as zero.
Unsupplied load cost is again assumed as 1,000 $/MW-month and the load for the five years is seen in Figure 9 with respective data presented in Table 9.
5.2. Sensitivity Analysis
To evaluate the proposed methodology, the mean operation cost for a single scenario using the proposed methodology was compared with the results for the same scenario using the Deterministic Dual Dynamic Programming (DDDP), considering 30 different initial reservoir storage levels.
In this first study 6, 11, and 21 discretizations, representing a variation of 20%, 10%, and 5% in the reservoir levels, were considered, respectively. Table 10 shows the values of the mean operation costs for these discretization levels for a five-year horizon, that is, 60 months.
Additionally, the mean operation cost of each discretization is compared with the mean cost by using the DDDP methodology, namely, $2,130,756.34. Table 3 shows the computational time required for each discretization level, using an Intel Core 2 Quad 2.5 GHz computer, with a 4 Gb RAM. From this result, 11 levels of discretization will be used to model the expected cost-to-go function in the next case study, as there is a slight difference to the DDDP cost, but with a superior computer performance than the 21 discretization solution.
5.3. Modeling the Expected Cost-to-Go Functions
Based on the latter study, this simulation considers 11 reservoir storage discretizations, in addition with 71 water inflow scenarios to model the expected cost-to-go function, in the backward phase. These scenarios are taken from the Brazilian hydrothermal system over 71 years of incremental water inflow recorded for each reservoir. Figure 10 shows a high, a low, and an average water inflow scenario, for the hydroplant 1.
Figure 11 shows the set of planes that models the expected cost-to-go function obtained from the proposed approach, in the last stage of the SDP problem (stage 60). Additionally, Figures 12 and 13 show the planes of stages 59 and 2, respectively, where there is a change in the morphology of these functions with the decreasing of the stage.
From the cost-to-go functions obtained in the backward phase, it is possible to calculate the expected total operation cost, in the forward phase, to different values of initial reservoir storage levels.
Although the mean operation cost takes into account all the scenarios, Table 11 shows the results for the three scenarios considered above: optimistic, pessimistic, and average.
The operation cost depends on the initial storage level, which is the amount of reservoir water storage in the first stage of the analysis. Three different initial storage levels were considered in this case study. The first is the upstream hydroplant with its reservoir almost full, about 96% of its maximum capacity of 34,116.00 hm3. The second and third cases show this same reservoir with an initial level around 58% and 29% of their maximum capacity. The last case is where the upstream reservoir is almost empty, near its minimum level of 5,447.00 hm3 as seen in Table 7. Additionally, as the downstream reservoir is smaller, this does not influence much in the operation cost.
In the optimistic scenario the initial reservoir storage level does not interfere in the operation cost. Due to the high incremental water inflow, the system can always store enough water to supply the load and recover the storage level, even when the system initially has a low reservoir storage level. Therefore, the three scenarios have the same expected operation costs, which is not the same for the other two scenarios.
For the average and pessimistic scenarios, in the first case, where the upstream reservoir (hydroplant 1) has a higher water storage, the expected cost to operate the system is significantly lower than the other two cases, with lower storage levels. This is as a result of the water level being lower, and it is necessary to turn on more expensive thermal generators to meet the load, impacting on the operation cost. This result reflects the influence of water inflow and reservoir storage in the long-term operation of hydrothermal systems.
5.4. Complete Simulation and Computational Platform
A Visual C++ platform was used to deal with the optimal hydrothermal system operation, using the SDP/Convex Hull approach. The last case study makes use of this platform.
The complete simulation includes 5 extra years in the backward phase to avoid the complete depletion of reservoirs in the last year of the given horizon. The same 11 discretization levels and 71-year records of water inflow were also considered. The initial storage levels were 13.3% and 31.8% for reservoirs 1 and 2, respectively.
The total operation cost, in the forward phase considering a single water inflow scenario, for the median of those 71-year records of water inflow, was expected to be $2,140,462.00 which is a high cost due to a large amount of additional thermal generation.
The problem took an approximate computing time of 8 minutes, with Figure 14 showing the amount of hydro- and thermal generation versus the total load for the 5 years horizon (2008 to 2012). The dark line represents the system demand, which is always supplied in this case study, as shown in the graph.
The advantage of optimizing the hydrothermal generation is shown in Figure 14, where the thermal generation level is analyzed. The optimization process uses an average thermal generation on a periodical basis, except for those months where the inflows are higher than expected. This agrees with the theoretical idea of committing the cheapest thermal units during the year, avoiding expensive thermal unit commitment, thus reducing the total cost.
Also, the computational system shows the hydroplants reservoir storage levels in Figures 15 and 16. In hydroplant 2, the reservoir is characterized by a rapid increase and decrease of its storage level. This is due to the fact that the hydroplant 1, which is located upstream, controls the river basin discharge and, also, this reservoir is significantly larger than the second one.
6. Conclusion
A linear piecewise approximation of expected cost-to-go functions of stochastic dynamic programming approach to the long-term hydrothermal operation planning using Convex Hull algorithms has been described in this study.
There is a tutorial system with a single hydroplant and a case study with two cascaded reservoir hydroplants and three thermal plants. The first study was a sensitivity analysis comparing the operation cost of the proposed methodology with the Deterministic Dual Dynamic Programming, in a single water inflow scenario. There is a convergence of both methodologies when an adequate number of discretizations are considered.
Another simulation, using the discretization parameter of the previous case study was carried out, to also show the acquisition of the expected cost-to-go functions. Finally, a complete simulation was described, for a complete long-term hydrothermal operation problem, for 2 hydroplants. These results show that it is possible to have the optimal operation of hydrothermal systems using the proposed methodology, when aggregated energy reservoirs are considered.
7. Future Work
Considering the results and the properties of the proposed methodologies, further ongoing studies are to extend it to a real complete system, like the Brazilian hydrothermal system with its four energy equivalent reservoir systems (South, Southeast/Center West, North, and Northeast). Besides, the transmission constraints should be included in order to have a more realistic approach.
To reduce the computational effort, an efficient sampling of the state space can be used, excluding states with equal costs or even states with reservoir storage levels that cannot be reached in a given stage. The use of parallel computing would also result in a significant reduction of this computational time.
Acknowledgments
The authors would like to thank the Brazilian National Council for Scientific and Technological Development—CNPq and FAPEMIG for financial support.