USE OF THE PAR(p) MODEL IN THE STOCHASTIC DUAL DYNAMIC PROGRAMMING OPTIMIZATION SCHEME USED IN THE OPERATION PLANNING OF THE BRAZILIAN HYDROPOWER SYSTEM

M. E. P. Maceira; J. M. Damázio

doi:10.1017/S0269964806060098

USE OF THE PAR(p) MODEL IN THE STOCHASTIC DUAL DYNAMIC PROGRAMMING OPTIMIZATION SCHEME USED IN THE OPERATION PLANNING OF THE BRAZILIAN HYDROPOWER SYSTEM

Published online by Cambridge University Press: 12 December 2005

M. E. P. Maceira and

J. M. Damázio

Show author details

M. E. P. Maceira: Affiliation:
Electric Power Research Center, CEPEL, Rio de Janeiro, Brazil, and, Rio de Janeiro State University, UERJ, Rio de Janeiro, Brazil, E-mail: elvira@cepel.br; damazio@cepel.br
J. M. Damázio: Affiliation:
Electric Power Research Center, CEPEL, Rio de Janeiro, Brazil, and, Rio de Janeiro State University, UERJ, Rio de Janeiro, Brazil, E-mail: elvira@cepel.br; damazio@cepel.br

Article contents

Abstract
1. INTRODUCTION
2. MODELING APPROACH
3. PROBLEM IDENTIFICATION
4. A PROPOSED SOLUTION
5. CONCLUSION
Acknowledgments
References

Rights & Permissions

Abstract

In September 2000, the Brazilian system dispatch and spot prices were calculated twice, using different inflow forecasts for that month, as in the last 5 days of August the inflows to the reservoirs in the South and Southeast regions changed 200%. The first run used a smaller forecasted energy inflow and the second used a higher energy inflow. Contrary to expectations, the spot price in the second run, with the higher energy inflow, was higher than the one found in the first run. This paper describes the problem, presents the special features of the PAR(p) model that allow the described behavior, and shows the solution taken to avoid the problem.

Information

Type: Research Article
Information: Probability in the Engineering and Informational Sciences , Volume 20 , Issue 1 , January 2006 , pp. 143 - 156

DOI: https://doi.org/10.1017/S0269964806060098 [Opens in a new window]
Copyright: © 2006 Cambridge University Press

1. INTRODUCTION

A chain of optimization models with different planning horizons and degrees of detail in system representation composes a computational system for the electric-energetic operation planning and programming of the Brazilian generating system [4]. The chain was proposed to the Brazilian Independent System Operator (ONS) and currently is being gradually implemented. ONS is responsible for central system optimization and dispatch according to clearly defined rules agreed upon by all industry members and approved by the regulatory body (ANEEL). The optimization models also compute the water values that form the basis for determining the Wholesale Energy Market (MAE) spot price.

On the top of the chain, the long-term operation planning model, called NEWAVE [3], applies a stochastic dual dynamic programming algorithm, which uses Bender's decomposition, to produce an operation strategy for each stage (month) of the planning period [6]. The operation strategy, given the system state at the beginning of the stage, produces generation targets for each of the four aggregated reservoirs, corresponding to the four Brazilian regions (South, Southeast, North, and Northeast, illustrated in Figure 1), and for the thermal plants, plus the energy interchanges among regions. The state system includes the stored energy in the beginning of each month for each of the four aggregated reservoirs or subsystems and information about “hydrologic trend” (i.e., the last p observed energy inflows in each of the four subsystems). This number, p, can vary from month to month and from subsystem to subsystem.

Brazilian subsystems.

In stochastic streamflow modeling, this corresponds to the use of the well-known and widely used periodic autoregressive PAR(p) model [9], which has been used as a forecast model for monthly flows and has been demonstrated in the comprehensive study of Noakes, McLeod, and Hippel [5], which considered several models.

Noakes et al. used 30 monthly unregulated river flow time series. The last 36 observations were omitted in the model identification and estimation phases of the study. One-step-ahead forecast statistics of the obtained forecasts in the omitted period were used to compare the models. The PAR(p) model obtained by identification of the model's order in each month using the partial autocorrelation function (PACF) and estimating the parameters using the Yule–Walker equations [9] was found to be the best model for the overall data in the study. Noakes et al. also commented, “Although other models may be more parsimonious, the PAR/PACF model gives the most adequate model due to the seasonal correlation effect” and “Besides being used for forecasting, PAR(p) models can also be employed for simulating hydrologic sequences” (i.e., to generate a great number of energy inflow scenarios) [2]. The good results obtained with the PAR/PACF modeling approach, both in one-step-ahead forecasts and simulation of drought sequences reasonably similar to the historical ones together with the easy mathematical tractability, recommended the use of this approach in the NEWAVE algorithm to model the energy inflow stochasticity. It was considered sufficient to use a maximum order p equal to 6, although the effective order varies for each month, according to the PAR/PACF modeling approach. This model is called PAR(6). The PAR(6) model generates energy inflow scenarios that will be used in the simulation process of the NEWAVE model and in the construction of the expected cost-to-go function of each month produced by the stochastic dual dynamic programming algorithm.

At the end of each month, the dispatch model is run by ONS, based on the inflow forecast for the next month. In addition to being used to determine the hydro and thermal generation, plus the energy interchanges among subsystems, the model calculates the spot price for each subsystem. The spot price is used for settling transactions in the Wholesale Energy Market and is therefore of great commercial importance.

In September 2000, the system dispatch and spot prices were calculated twice, using different inflow forecasts for that month, as in the last 5 days of August the inflows to the reservoirs in the South and Southeast regions changed 200%. The first run used a smaller forecasted energy inflow and the second used a higher energy inflow. Contrary to expectations, the spot price in the second run, with the higher energy inflow, was higher than the one found in the first run.

The analysis of the case made by [7] suggested that the problem occurred because of an error in the calculation of Bender's cut that constitutes the expected cost-to-go-function of each month.

The problem was further investigated and it was detected that the expected energy inflows produced by the PAR(p) model in February 2001 decrease sharply when the energy inflow in September 2000 was increased and this was the reason of the increase in the spot price. The suggestion in [8] proved not be the case.

The objective of this work is to describe the problem that occurred in September 2000 in the operation planning of the Brazilian hydropower system, present the special features of the PAR(p) model that allow the described behavior, and show the solution taken to avoid the problem.

2. MODELING APPROACH

The operation strategy minimizes the expected value of the operation cost during the planning period, which is composed of fuel costs plus penalties for failure in load supply. If the inflow volumes are known at the beginning of each stage, stochastic dual dynamic programming [2,6,7], represented by the following backward recursive equation, can solve the operation dispatch problem:

subject to the following:

Storage balance equation in each subsystem aggregated reservoir

Demand supply equation in each subsystem k

Bounds in storage in each subsystem aggregated reservoir

Maximum hydro production in each subsystem aggregated reservoir

Lower bounds on total outflow in each subsystem aggregated reservoir

Maximum generation in each thermal plant

Flow limits among subsystems

Set of multivariate linear constraints representing the cost-to-go function (Bender's cut)

where

The periodic autoregressive model, denoted by PAR(p), can be written as

where Z_t,t = 1,2,…, is the seasonal time series with period s (in the case of the monthly time series, s is equal to 12). The time index t may be regarded as a function of the year T (T = 1,…,N) and the season m (m = 1,…,s): t = (T − 1) s + m. μ_m is the mean of period of m, σ_m is the standard deviation of period m, and a_t is the time uncorrelated series, which is independent of Z_t, and it has also zero mean and variance equal to σ_a^2(m).

Take π_t, π_d,t, and λ_t as the Lagrange multipliers associated to the storage balance equation, demand supply equation, and cost-to-go multivariate linear constraints, respectively, of the dispatch problem of stage t. After the solution of each operation dispatch problem, the Bender's cut coefficients π_v,t,π_A1,t,…,π_Ap,t that will be added to the expected cost-to-go function of stage (t − 1), are calculated deriving the operation cost in relation to the energy stored in the beginning of stage t and the energy inflows in stages (t − 1),…,(t − p), obtaining

3. PROBLEM IDENTIFICATION

When, in September 2000, the system dispatch and spot prices were calculated twice, using different inflow forecasts for September and the consolidated and known energy inflows for August, July, and so on, the spot price associated with a lower forecasted energy inflow could be smaller than that associated with the higher forecasted energy inflow only if the π_A1,t Bender's cut coefficient is positive.

Next, the Bender's cut coefficient calculation responsible for a π_A1,t positive in the dispatch operation problem of September 2000 will be detailed. The adopted horizon of the backward recursive simulation was 10 years, starting in December 2009 and ending in September 2000. To illustrate the problem, it will be sufficient to carry out the application of the recursive equation (1) from February 2001 to September 2000 considering only one energy inflow scenario. Table 1 shows the PAR(p) coefficients of the Southeast subsystem from February to October.

PAR(p) Coefficients of the Southeast Subsystem

Table 2 shows the simplex multipliers associated with storage balance, demand supply equations, and the active linear constraint representing the cost-to-go function, in addition to the Bender's cut coefficients of the Southeast subsystem from February 2001 to October 2000.

Operation Dispatch Parameters

With the values of Table 2a, it is possible to calculate the Bender's cut coefficients that will be aggregated to the expected cost-to-go function of January 2001:

With the values of Table 2b, it is possible to calculate the Bender's cut coefficients that will be aggregated to the expected cost-to-go function of December 2000:

With the values of Table 2c, it is possible to calculate the Bender's cut coefficients that will be aggregated to the expected cost-to-go function of November 2000:

With the values of Table 2d, it is possible to calculate the Bender's cut coefficients that will be aggregated to the expected cost-to-go function of October 2000:

With the values of Table 2e, it is possible to calculate the Bender's cut coefficients that will be aggregated to the expected cost-to-go function of September 2000:

Now, the Bender's cut that will be aggregated to expected cost-to-go function of September 2000 can be written as

It can be observed that the Bender's cut coefficient associated to the forecasted energy inflow A_t is, in this case, positive. The appearance of the positive coefficient in the Bender's cut is due to the expression of the expected energy inflow of February 2001 for the Southeast subsystem as a function of the state at September 2000 and includes the inflow at September 2000 times a negative constant. This can be shown by back-substitution of the equations of the January inflow, December inflow, November inflow, October inflow, and September inflow in the equation of the February inflow, as is shown in the following.

Considering X_t as (Z_t − μ_m), the regression of February in the PAR(p) model can be written as

Substituting the known regression of January in the PAR(p) model (X_Jan = 0.934X_Dec) in the above expression, it can be written as

Regrouping common terms, the regression of February can be written as

In the same way, the December regression is given by X_Dec = 0.971 X_Nov. Substituting the following expression in the previous one, we obtain

Substituting the known PAR(p) model regression of November (X_Nov = 0.816X_Oct) in the above expression results in

Once more, the October regression is given by X_Oct = 0.470X_Sept − 0.151X_Aug + 0.581X_July. Substituting this in the previous expression and regrouping the terms, the regression of February 2001 inflow as a function of the September inflow forecast and past inflows is obtained:

This regression shows a negative coefficient multiplying the September energy inflow forecast. It can be noted that the negative September coefficient, which showed up after these back-substitutions, has its roots in negative parameters in the February model for the months 2-lagged (December), 4-lagged (October), and 5-lagged (September).

4. A PROPOSED SOLUTION

The identification procedure for the February model was reappraised. The PACF of this month (Fig. 2) shows a significant value at lag 1, nonsignificant values from lag 2 to lag 5, and a significant value at lag 6. Using the adopted classical identification rule, “choose p as the largest significant order i so that the partial correlations for k > i are not significant,” this last significant value yielded a model of order 6. On the order hand, even a pure Box and Jenkins's identification rule application [1] would enhance the classical identification rule, by neglecting not very strong significant partial correlations. This would turn the model more parsimonious without increasing its residual sum of squares very much. After a general reappraisal of the identifications for the Brazilian subsystems, it was discovered that the same pattern occurred in two other months of the South subsystem and in one month for the North and Northeast subsystems.

February PACF.

Taking out the occurrence of runs of nonsignificant partial correlations before a significant one, as in these particular cases, in general one may state that, with the exception in the cases with flows affected by snowmelt, there is no physical reason for the inflow in any month to affect any following inflows by a negative coefficient. This would pose a positive constraint in all of the coefficients of PAR(p) models. On the other hand, with a pure statistical point of view, the obtained negative coefficients are the best representation of the correlation pattern revealed by the data, and the literature on PAR(p) modeling of streamflows is full of examples of estimated negative coefficients; examples can be found in [5] and [10]. As a compromise solution, it was proposed to maintain the classical rule of identification but to reduce the order of the PAR(p) model of each stage whenever the first coefficient resulting from any of 12 back-substitutions of the equation's model turns out to be negative. This model was called PAR(6a). The result was equivalent to the above-cited approach of considering those, maybe spurious, significant partial correlations, as insignificant and considering the more parsimonious modeling. It was shown that this PAR(6a) model can also reproduce the historical occurrence of severe droughts. Table 3 compares the order selected from models PAR(6) and PAR(6a).

Comparison of Order Selected from Models PAR(6) to PAR(6a), Number of Autoregressive Coefficients

Figure 3 compares the historical (70 years) and synthetic (2000 years) average energy inflow produced by the PAR(6) and PAR(6a) models for the Southeast subsystem starting from different September energy inflows. Two sets of 2000 synthetic energy inflows are presented. The solid curve was obtained using the near-average energy inflow for September of the first run dispatch operation model in September 2000. The long-dashed curve was obtained using the higher energy inflow for September of the second run dispatch operation model in September 2000. It can be seen that the parsimonious model, PAR(6a), presents a more reliable behavior.

Historical × synthetic average energy inflows of the Southeast models: (a) PAR(6) and (b) PAR(6a).

5. CONCLUSION

The increased spot price with higher-energy inflow forecast that occurred in September 2000 in the Brazilian dispatch operation mode is not due to errors in formulation or implementation in the stochastic dual dynamic programming algorithm. Otherwise, it has been shown that the use of the classical approach of the PAR(p) model can potentially produce this problem.

The proposed solution maintains the classical rule of identification but reduces the order of the model whenever the first coefficient resulting from 12 back-substitutions of the regressions turns out to be negative. These more parsimonious models presented a more reliable behavior even when applied to the September 2000 Brazilian dispatch problem.

Acknowledgments

The authors would like to thank L.A. Terry, F.S. Costa, and A.C.G. Melo for their very important contributions to the development of this work.

References

REFERENCES

Box, G.E.P. & Jenkins, G.M. (1970). Time Series Analysis: Forecasting and Control. Littleton, CO: Holden-Day.

Maceira, M.E.P. & Bezerra, C.V. (1997). Stochastic streamflow model for hydroelectric systems. In 5th Probabilistic Methods Applied to Power Systems—PMAPS-1997.

Maceira, M.E.P., Mercio, C.B., Gorenstin, B.G., Cunha, S.H.F., Suanno, C., Sacramento, M.C., & Kligerman, A.S. (1998). Energy evaluation of the North/Northeastern and South/Southeastern interconnection with NEWAVE model. In VI Symposium of Specialists in Electric Operational and Expansion Planning—SEPOPE.

Maceira, M.E.P., Terry, L.A., Damazio, J.M., Costa, F.S., & Melo, A.C.G. (2002). Chain of optimization models for setting the energy dispatch and spot price in the Brazilian system. In Power System Computation Conference—PSCC'02.

Noakes, D.J., McLeod, A.I., & Hipel, K.W. (1985). Forecasting monthly riverflow times series. International Journal of Forecasting 1: 179–190.Google Scholar

Pereira, M.V.F. (1989). Optimal stochastic operations scheduling of large hydroelectric systems. International Journal of Electric Power and Energy Systems 11(3): 161–169.Google Scholar

Pereira, M.V.F. & Pinto, L.M.V.G. (1985). Stochastic optimization of a multireservoir hydroelectric system: A decomposition approach. Water Resources Research 21(6): 779–792.Google Scholar

PSR Consultoria. (2000). Analysis of spot price projection for September and October of 2000. Technical report, Draft 2. PSR Consultoria, Rio de Janeiro.

Salas, J.D., Delleur, J.W., Yevjevich, V., & Lane, W.L. (1980). Applied Modeling of Hydroeletric Series. Water Resources Publications, San Francisco.

Thompstone, R.M., Hipel, K.W., & Mcleod, A.I. (1984). Forecasting quarter-monthly riverflow. Technical report TR-84-07, Department of Statistical & Actuarial Sciences, University of Western Ontario, London, Ontario, Canada.