Transcript talk - University of Auckland

Stochastic processes for hydrological
optimization problems
Geoffrey Pritchard
University of Auckland
Prologue: Iterated Function Systems
Consider the following Markov process in the plane:
  0.2  0.4  X t   0.5 
     , 50% chance
 
0.7 0.6  Yt   0.6 


X
 t 1  

  
 Yt 1    0.8 0.5  X   0.4 
 t     , 50% chance
 
  0.1  0.9  Yt   0.7 
(Each step is randomly chosen from a finite list of affine
transformations, independently of previous steps.)
(Is this useful for anything, besides (maybe) computer graphics?)
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.
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.
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
Time/information structure
Week t-1
E
Week t
Week t+1
min (present cost) + E[ future cost ]
s.t. satisfy demand, etc. with
stored energy + random inflow Xt
• 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.
Optimization
Week t-1
gt-1(y) =
E
Week t
Week t+1
minu (present cost)(u) + gt(u)
s.t. satisfy demand, etc. with
(stored energy)(y) + random inflow Xt
• Let gt(u) = expected cost of consequences after week t of doing u in week t.
• Essential observation is that g t  g t 1 is convexity-preserving.
– so all optimizations can be of convex functions, i.e. tractable.
• Computationally, convex subproblems -> linear programs.
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
Inflow data
Waitaki catchment above Benmore dam, weekly, 1948-2010
Strong seasonal dependence
– 3:1 ratio between midsummer high/midwinter low.
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).
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!
De-seasonalization
inflow
Xt
Qt 
gt
via regression:
log X t  log gt   t
A convenient
normalization, but does
not make (Qt) stationary!
Suggestion: an autoregressive model
The AR(1) model
ACF of (Qt)
log Qt  t log Qt 1   t
that is,
Qt  Qtt1 exp(  t )
seems reasonable.
(Life should be so simple.)
Stagewise independence
Week t-1
gt-1(y) =
E
Week t
Week t+1
minu (present cost)(u) + gt(u)
s.t. satisfy demand, etc. with
(stored energy)(y) + random inflow Xt
• 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.
From independent to Markov inflows
Week t
Week t-1
gt-1(y) =
E
Week t+1
minu (present cost)(u) + gt(u)
s.t. satisfy demand, etc. with
(stored energy)(y) + (inflow)(y, Wt)
• 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.
LP-compatible autoregressive processes
• We’re allowed a process with
Qt  H t (Qt 1 ,Wt )
with (Wt ) an independent sequence, and H t (  ,Wt ) a linear function.
• But what we had in mind was
Qt  Qtt1 exp(  t )
which is nonlinear.
• It’s concave, though, so admits a
piecewise linear approximation.
LP-compatible autoregressive processes
• We’re allowed a process with
Qt  H t (Qt 1 ,Wt )
with (Wt ) an independent sequence, and H t (  ,Wt ) a linear function.
• But what we had in mind was
Qt  Qtt1 exp(  t )
which is nonlinear.
• It’s concave, though, so admits a
piecewise linear approximation.
LP-compatible autoregressive processes
• We’re allowed a process with
Qt  H t (Qt 1 ,Wt )
with (Wt ) an independent sequence, and H t (  ,Wt ) a linear function.
• But what we had in mind was
Qt  Qtt1 exp(  t )
which is nonlinear.
• It’s concave, though, so admits a
piecewise linear approximation.
One linear piece
Approximate the model
Qt  Qtt1 exp(  t )
by linearizing
q  q t
about
q 1
Qt  (1  t (Qt 1  1)) exp(  t )
How to fit this?
Inference
1. Auto-regression on log-inflows:
log Qt  t log Qt 1   t
(ignores the linear approximation step)
2. Auto-regression on Qt-1 :
Qt  1  t (Qt 1  1)   t
(ignores the structure of errors)
3. Or we could do it right: a max-likelihood fit on the actual model:
Qt  (1  t (Qt 1  1)) exp(  t )
Stochastic dual dynamic programming (SDDP)
Week t
Week t-1
gt-1(y) =
E
Week t+1
minu (present cost)(u) + gt(u)
s.t. satisfy demand, etc. with
(stored energy)(y) + random inflow Xt
• 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.
The importance of being discrete
Week t-1
gt-1(y) =
Ss
Week t
Week t+1
minu (present cost)(u) + gt(u)
s.t. satisfy demand, etc. with
(stored energy)(y) + random inflow Xt
ps
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)
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.
A catalogue of errors
Our efforts to model inflows have incurred
• model mis-specification error
– inflows might not really be AR-1
discrete (data)
• inferential sampling error (finite data)
– parameters may be wrong
continuous (AR-1 model)
• sample average approximation error
– optimization is vs. a discrete
approximation of inflow process
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?
The final form required
Model for inflow Qt in week t :
Qt  (1  t (Qt 1  1)) exp(  t )
(t discretely distributed)
Or more generally
Qt  Rt  St Qt 1
- where (Rt, St) is chosen at random from a small collection of
(seasonally-varying) scenarios.
A linear iterated function system (IFS)
Markov process.
The final form required
Model for inflow Qt in week t :
Qt  (1  t (Qt 1  1)) exp(  t )
(t discretely distributed)
Or more generally
Qt  Rt  St Qt 1
- where (Rt, St) is chosen at random from a small collection of
(seasonally-varying) scenarios.
A linear iterated function system (IFS)
Markov process.
Linear IFS Markov inflow model
Fitting a model to data: quantile regression
• Have data xi and yi for i=1,…n
y
x
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
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
• Quantile regression:
choose scenario sk() to
minimize Si rk( yi – sk(xi) )
for a suitable loss function rk().
x
Quantile regression fitting
• For a scenario at quantile t (0 < t < 1) ,
rt is the loss function
t1
t
• For each scenario, the quantile regression problem is a linear program.
Fitting a model to data: quantile autoregression
Quantiles (t)
Scenario probabilities
0.02
0.1
0.2
0.35
0.5
0.06
0.07
0.125
0.15
0.15
0.65
0.15
0.8
0.125
0.9
0.98
0.07
0.06
For each of a fixed collection of quantiles, fit a scenario (Rt, St) by
quantile regression:
X t  Rt  St X t 1
Scenario probabilities determined from quantiles; can be unequal.
Fitting a model to data: quantile autoregression
X t  Rt  St X t 1
• 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.
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,

CRPS( F , y)   ( F ( x)  1x< y ) dx
2

F
y
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
y
• Choose sk() and pk to
minimize


j CRPS k pk sk ( x j ) , y j 
x
• This is the most computationally challenging method
(global optimization, not LP or least-squares).
Fitting a model to data: CRPS M-estimation
X t  Rt  St X t 1
• Scenarios may cross.
• Scenario probabilities are
optimized in model fitting, instead
of being arbitrarily chosen
• Quite small numbers of scenarios
seem possible.
Multivariate inflow models
 North Island inflow 
  Rt  St X t 1
X t  
 South Island inflow 
• 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.
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).
Results – optimal strategy
Lost load
Spill
(MW, probability) (MW, probability)
Energy
price
($/MWh)
Inflow model
No. scenarios
per stage
quantile
regression
16
9.4, 28%
2.9,
6%
251
autoregressive,
resampled errors
63
13.3, 37%
17.1, 15%
296
autoregressive,
lognormal errors
63
6.8, 23%
6.0, 9%
207
independent
(uncorrected)
16
1.59,
9%
0.14,
1%
(Quantities are expected averages over full time horizon;
probabilities are for any shortage/spill within time horizon.)
112