Generating time series for electricity demand

Question

I'm working on a project, trying to model electricity consumption and generation from a bunch of PV generators to see how much of the demand can be satisfied by the production in the "town". Most models use a fixed demand and production profile in an hourly time-series.

I'd like to see how variations in both production and demand would influence the outcome. So basically a Monte Carlo simulation of the whole system.

Is there a relatively easy way to vary a given demand and production profile with a sort of random factor? Obviously the diurnal and annual variations should still be similar to the "measured" dataset.

I played around with something called Iterative Amplitude Adapted Fourier Transform, but the variations are totally out of whack (maybe because I didn't adjust the parameters correctly). Sadly my education in statistics is not very deep.

Thanks for any suggestion!

Answer 1

The easiest way would probably to take an electricity demand (or production) time series, fit an appropriate model, randomly permute the residuals and add these back to the fitted values. You can inflate the permuted residuals to make the simulated series less similar to the original series.

Alternatively, don't resample the residuals, but instead draw from an appropriate parametric distribution. However, given the amount of data you will likely be working with, I don't this this will yield very different results from the permutation approach.

Here is an example. A good model to fit electricity data is tbats, which can model multiple-seasonalities. The top left panel is the original series; the rest are resampled simulations. Note how the intra-daily and intra-weekly patterns look similar, but the simulated series are more noisy. You can tweak the scaling-factor to govern this.

library(forecast)
taylor.fit <- tbats(taylor)

scaling.factor <- 100000
opar <- par(mfrow=c(3,3),mai=c(.4,.4,0,0))
    plot(taylor,xlab="",ylab="",ylim=c(15000,44000))
    for ( ii in 1:8 ) {
        set.seed(ii)
        plot(fitted(taylor.fit)+scaling.factor*sample(residuals(taylor.fit),length(taylor)),
            xlab="",ylab="",ylim=c(15000,44000))
    }
par(opar)