Modéliser la consommation électrique

A partir de nos prévisions pour la température, on peut tenter de prévoir la consommation électrique. Rappelons que la série de consommation électrique ressemble à ca

plot(electricite[passe,"Load"],type="l")

On peut tenter un modèle assez simple, où la consommation à la date Y_t est fonction d’une tendance linéaire a+bt, de la position dans l’année (sous une forme non-paramétrique), de la température T_t (sous une forme non-linéaire, disons un polynôme de degré 3), de la températeure en t-1, T_{t-1} et de l’indice de production industriel IPI_t (sous une forme linéaire pour faire simple).

model1 = lm(Load~1+Time+as.factor(NumWeek)+
poly(Temp,3)+Temp1+IPI,data=electricite[passe,])

L’idée de la fonction polynomial pour l’impact de la température venait du graphique suivant (où on représente la série de la consommation, une fois enlevée la tendance linéaire)

On peut aussi supposer une forme autorégressive, où Y_{t} sera fonction de Y_{t-1}

model2 = lm(Load~1+Load1+Time+as.factor(NumWeek)+
poly(Temp,3)+Temp1+IPI,data=electricite[passe,])

On peut alors tenter de faire de la prévision. Le soucis avec ce second modèle, c’est la composante autorégressive. Pour prévoir Y_{t+h}, il faut utiliser la prévision faite en t+h-1, \hat Y_{t+h-1}.

base_prevision=data.frame(
Load_prev = NA,  Load1=c(electricite[length(passe),"Load"],rep(NA,110),
Time = electricite[futur,"Time"],
NumWeek = electricite[futur,"NumWeek"],
Temp = MODELES[,"y4"], 
Temp1 = c(electricite[length(passe),"Temp"],
MODELES[1:110,"y4"]),
IPI = electricite[futur,"IPI"])
for(t in 1:110){
 p = predict(model2,newdata=base_prevision)[t]
 base_prevision[t,"Load_prev"] = p
 base_prevision[t+1,"Load1"] = p}

On peut alors comparer la prévision \hat Y_{t} avec l’observation Y_{t}.

plot(1:111,electricite[futur,"Load"],type="l")
lines(1:111,base_prevision[,"Load_prev"],col="red")

On est bon les mois d’été (on capture bien le pic de la première quinzaine du mois d’août) par contre, on sous-estime fortement la consommation l’hiver. Et c’est gênant…

Pour finir, on peut oublier le passage par les variables explicatives, et tenter directement un modèle de série temporelles.

plot(electricite[passe,"Load"],type="l")

Le soucis, c’est le côté hétéroscédastique de la série, avec une pente sur le minimum plus faible que sur le maximum.

n=length(passe)
plot(1:n,electricite[passe,"Load"],type="l")
m=aggregate(electricite[passe,"Load"],
by=list(as.factor(electricite[passe,"Year"])),min)
M=aggregate(electricite[passe,"Load"],
by=list(as.factor(electricite[passe,"Year"])),max)
vm=which(electricite[passe,"Load"]%in%m$x)
vM=which(electricite[passe,"Load"]%in%M$x)
points(sort((1:n)[vm]),m$x,pch=19,col="blue")
points(sort((1:n)[vM]),M$x,pch=19,col="blue")
xm=((1:n)[vm])
xM=((1:n)[vM])
regm=lm(m$x~xm,col="blue")
regM=lm(M$x~xM,col="blue")
abline(regm,lty=2)
abline(regM,lty=2)

Une solution classique (et simple) est de passe au logarithme

plot(electricite[passe,"Load"],type="l",log="y")
y=electricite[passe,"Load"]
z = log(y)
plot(z,type="l")
B = data.frame(z=z,t=1:length(z))
abline(lm(z~t,data=B),lty=2,col="red")

Il faut ensuite enlever la tendance linéaire, pour stationnariser la série

z = residuals(lm(z~t,data=B))
Z = ts(z,start=1996,frequency=52)
modelz1 = arima(Z,order = c(0,0,0), 
seasonal = list(order = c(1,0,0)))
modelz1 = arima(Z,order = c(4,0,0), 
seasonal = list(order = c(1,0,1)))

Ce premier modèle est stationnaire, sans racine unité. On peut tenter d’introduire une racine unité saisonnière

modelz2 = arima(Z,order = c(0,0,0), 
seasonal = list(order = c(0,1,0)))
modelz2 = arima(Z,order = c(1,0,3), 
seasonal = list(order = c(1,1,0)))

Enfin, le dernier sera un peu plus simple

modelz3 = arima(Z,order = c(1,0,0), 
seasonal = list(order = c(2,0,0)))

On stocke alors toutes les prévisions dans une base

DONNEES_u = data.frame(
  u1 = predict(modelz1,n.ahead = 111)$pred,
  u2 = predict(modelz2,n.ahead = 111)$pred,
  u3 = predict(modelz3,n.ahead = 111)$pred)

On rajoute alors la tendance linéaire à la prévision du résidu

reg = lm(z~t,data=B)
dfutur = data.frame(t=futur)
DONNEES_z = data.frame(
  z1 = DONNEES_u$u1+predict(reg,newdata=dfutur),
  z2 = DONNEES_u$u2+predict(reg,newdata=dfutur),
  z3 = DONNEES_u$u3+predict(reg,newdata=dfutur)

On a besoin de la variance pour la prévision. En effet, classiquement, la prévision sera l’espérance de notre grandeur d’intérêt. Ici, on a fait un modèle linéaire sur \log Y, i.e. \log Y \sim\mathcal{N}(\mu,\sigma^2) et donc \mathbb{E}[Y]=\exp\left(\mu+\sigma^2/2\right)

sigma=summary(reg)$sigma
DONNEES_seu = data.frame(
  seu1=sqrt(predict(modelz1,n.ahead = 111)$se^2+sigma^2),
  seu2=sqrt(predict(modelz2,n.ahead = 111)$se^2+sigma^2),
  seu3=sqrt(predict(modelz3,n.ahead = 111)$se^2+sigma^2))

On fait ici l’hypothèse que les estimateurs des prévisions des deux modèles (la tendance linéaire et celle du modèle autorégressif) sont indépendants et donc on peut sommer les termes de variance. Aussi, les prévisions pour Y sont

DONNEES_y = data.frame(
  y1 = exp(DONNEES_z$z1+1/2*DONNEES_seu$seu1^2),
  y2 = exp(DONNEES_z$z2+1/2*DONNEES_seu$seu2^2),
  y3 = exp(DONNEES_z$z3+1/2*DONNEES_seu$seu3^2))

On peut représenter les prévisions des trois modèles (avec les observations)

plot(futur,electricite[futur,"Load"],type="l")
lines(futur,DONNEES_y$y1,col="red")
lines(futur,DONNEES_y$y2,col="blue")
lines(futur,DONNEES_y$y3,col="purple")

On peut comparer avec la prévision précédente,

lines(futur,base_prevision[,"Load_prev"],
col="orange")

On y perd un peu l’été, par contre, on s’est un peu amélioré l’hiver.


One thought on “Modéliser la consommation électrique”

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.