Stochastic processes for hydrological optimization problems

1. Stochastic processes for hydrological optimization problems Geoffrey Pritchard University of Auckland

2. Prologue: Iterated Function Systems Consider the following Markov process in the plane: (Each step is randomly chosen from a finite list of affine transformations, independently of previous steps.)

13. (Is this useful for anything, besides (maybe) computer graphics?)

14. Hydro scheduling – an optimal control problem random inflows state variables: stored energy control variables: outflows Control water releases over time to maximize value. • As a workably competitive market would do.

15. Hydro scheduling – an optimal control problem statistics OR How not to do it: 1. Develop a stochastic model of inflows. 2. Optimize releases versus the given inflow-generating process.

16. Why develop a model of inflows? Why not just use the historical data non-parametrically? • Small dataset. e.g. autumn 2014 : - Mar ~ 1620 MW - Apr ~ 2280 MW - May ~ 4010 MW Past years (if any) with this exact sequence are not a reliable forecast for June 2014. • A model allows events to be more extreme than anything in the data The worst event ever observed is not the worst possible

17. Time/information structure Week t Week t+1 Week t-1 min (present cost) + E[ future cost ] s.t. satisfy demand, etc. with stored energy + random inflow Xt E • Each stage subproblem is a random optimization problem. • Stage t subproblem is solved with knowledge of Xt, but not of the future. • Weekly stages might be good, continuous time worse.

18. Optimization Week t Week t+1 Week t-1 minu (present cost)(u) + gt(u) s.t. satisfy demand, etc. with (stored energy)(y) + random inflow Xt E gt-1(y) = • Let gt(u) = expected cost of consequences after week t of doing u in week t. • Essential observation is that is convexity-preserving. • so all optimizations can be of convex functions, i.e. tractable. • Computationally, convex subproblems -> linear programs.

19. Model me this... (the upper Waitaki catchment) • ~ 20% of NZ electricity derives from precipitation in this region • Rainfall + summer snowmelt from Southern Alps. Tekapo A Tekapo B Ohau A Ohau B Ohau C Benmore Aviemore Waitaki

20. Inflow data Waitaki catchment above Benmore dam, weekly, 1948-2010 Strong seasonal dependence • 3:1 ratio between midsummer high/midwinter low.

21. Serial dependence • Weather patterns persist • increases probability of shortage/spill. • Typical correlation length ~ several weeks (but varying seasonally). • convenient for optimization (cf. e.g. Brazil).

22. Extreme values Hydro-scheduling is sensitive to extremes of inflow (in both tails). • Low inflow -> reservoirs run dry (the most momentous thing that can happen) • High inflow -> economic loss (spill); removes risk of shortage. • Beware discrete approximations to the distribution!

23. De-seasonalization inflow via regression: A convenient normalization, but does not make (Qt) stationary!

24. Suggestion: an autoregressive model The AR(1) model ACF of (Qt) that is, seems reasonable. (Life should be so simple.)

25. Stagewise independence Week t Week t+1 Week t-1 minu (present cost)(u) + gt(u) s.t. satisfy demand, etc. with (stored energy)(y) + random inflow Xt E gt-1(y) = • Stage t subproblem is solved with knowledge of Xt, but not of the future. • stagewise independence, i.e. (Xt) an independent sequence. • Inflows are not stagewise independent. • Suggested model is Markov.

26. From independent to Markov inflows Week t Week t+1 Week t-1 minu (present cost)(u) + gt(u) s.t. satisfy demand, etc. with (stored energy)(y) + (inflow)(y, Wt) E gt-1(y) = • Make inflow a function of • what happened last week (y) • a random innovation Wt – with (Wt ) independent That’ll work – if we can express it as a linear program.

27. LP-compatible autoregressive processes • We’re allowed a process with with an independent sequence, and a linear function. • But what we had in mind was which is nonlinear. • It’s concave, though, so admits a piecewise linear approximation.

28. LP-compatible autoregressive processes • We’re allowed a process with with an independent sequence, and a linear function. • But what we had in mind was which is nonlinear. • It’s concave, though, so admits a piecewise linear approximation.

29. LP-compatible autoregressive processes • We’re allowed a process with with an independent sequence, and a linear function. • But what we had in mind was which is nonlinear. • It’s concave, though, so admits a piecewise linear approximation.

30. One linear piece Approximate the model by linearizing about How to fit this?

31. Inference 1. Auto-regression on log-inflows: (ignores the linear approximation step) 2. Auto-regression on Qt-1 : (ignores the structure of errors) 3. Or we could do it right: a max-likelihood fit on the actual model:

33. Stochastic dual dynamic programming (SDDP) Week t Week t+1 Week t-1 minu (present cost)(u) + gt(u) s.t. satisfy demand, etc. with (stored energy)(y) + random inflow Xt E gt-1(y) = • The leading algorithm for problems of this type. • Essential step (backward pass): • evaluate the expectation for given y, using current estimate of gt. • Use dual variables from optimization to form a cut (linear lower bound), which improves estimate of gt-1.

34. The importance of being discrete Week t Week t+1 Week t-1 minu (present cost)(u) + gt(u) s.t. satisfy demand, etc. with (stored energy)(y) + random inflow Xt Ss ps gt-1(y) = s • If random elements have a discrete joint distribution: • solve the optimization problem for each atom. • Otherwise, need a (Monte Carlo?) discrete approximation • with not too many discrete scenarios, please (computation time is (at least) proportional to number of scenarios)

35. Sample average approximation (SAA) Objective function is an expectation, over a continuous distribution. Only way to evaluate it is by Monte Carlo sampling. Fix a sample, optimize the resulting approximation.

36. A catalogue of errors Our efforts to model inflows have incurred • model mis-specification error • inflows might not really be AR-1 • inferential sampling error (finite data) • parameters may be wrong • sample average approximation error • optimization is vs. a discrete approximation of inflow process discrete (data) continuous (AR-1 model) discrete (SAA approx to transition kernel) Can we obtain, in one step, a good representation of the data by a model of the final form required?

37. The final form required Model for inflow Qt in week t : (et discretely distributed) Or more generally - where (Rt, St) is chosen at random from a small collection of (seasonally-varying) scenarios. A linear iterated function system (IFS) Markov process.

38. The final form required Model for inflow Qt in week t : (et discretely distributed) Or more generally - where (Rt, St) is chosen at random from a small collection of (seasonally-varying) scenarios. A linear iterated function system (IFS) Markov process.

39. Linear IFS Markov inflow model

40. Fitting a model to data: quantile regression • Have data xi and yi for i=1,…n y x

41. Fitting a model to data: quantile regression • Have data xi and yi for i=1,…n • Want to represent the distribution of y|x by finitely many scenarios. y x

42. Fitting a model to data: quantile regression • Have data xi and yi for i=1,…n • Want to represent the distribution of y|x by finitely many scenarios. • Quantile regression: choose scenario sk() to minimize Sirk( yi – sk(xi) ) for a suitable loss function rk(). y x

43. Quantile regression fitting • For a scenario at quantilet (0 < t < 1) , rt is the loss function t-1 t • For each scenario, the quantile regression problem is a linear program.

44. Fitting a model to data: quantile autoregression Quantiles (t) 0.02 0.1 0.2 0.35 0.5 0.65 0.8 0.9 0.98 Scenario probabilities 0.06 0.07 0.125 0.15 0.15 0.15 0.125 0.07 0.06 For each of a fixed collection of quantiles, fit a scenario (Rt, St) by quantile regression: Scenario probabilities determined from quantiles; can be unequal.

45. Fitting a model to data: quantile autoregression • Scenarios should not cross. • Dependence (slope St) can vary across the probability distribution. High-flow scenarios differ in intercept (current rainfall). Low-flow scenarios differ mainly in slope. Extreme scenarios have their own dependence structure.

46. Continuous ranked probability score (CRPS) • A method of judging the merit of a prediction made in the form of a probability distribution. • Given prediction distribution F and actual outcome y, F y

47. Fitting a model to data: CRPS M-estimation • CRPS can also be used as an estimation method for multi-scenario regression. • Given x, scenarios for y are s1(x) … sm(x) with probabilities p1 ... pm • Choose sk() and pk to y x • This is the most computationally challenging method (global optimization, not LP or least-squares).

48. Fitting a model to data: CRPS M-estimation • Scenarios may cross. • Scenario probabilities are optimized in model fitting, instead of being arbitrarily chosen • Quite small numbers of scenarios seem possible.

49. Multivariate inflow models • Need to capture spatial as well as temporal correlations. • Generalize models: • Autoregressive: need discrete approx. to multivariate error. • Quantile regression: no natural generalization of quantile. • CRPS M-estimation: generalization to energy score.

50. A test problem • Challenging fictional system based on Waitaki catchment inflows. • Storage capacity 1000 GWh (cf. real Waitaki lakes 2800 GWh) • Generation capacity 1749 MW hydro, 900 MW thermal • Demand 1550 MW, constant • Thermal fuel $50 / MWh, VOLL $1000 / MWh • Test problem: a dry winter. • 35 weeks (2 April – 2 December) • Initial storage 336 GWh • Initial inflow 500 MW (~56% of average) • Solved with Doasa 2.0 (EPOC’s SDDP code).