Tag Archives: électricité

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.

Devoir, Méthodes de Prévisions (consommation électrique)

Pour le projet de l’UE de prévision, je demande une prévision de la consommation d’électricité. Pour construire les modèles, je donne accès à un gros fichier qui donne, pour toute la France, la consommation d’électricité, entre 2012 et 2014, heure par heure.

  • comme en cours, la première étape sera de proposer une série de 52 prévisions pour chaque semaine de 2015,
  • la seconde étape sera de proposer deux séries de 31 prévisions, pour chaque jour d’août 2015 et de décembre 2015
  • enfin, la dernière étape sera une prévision de 168 prévisions, pour chaque heure de la semaine du 19 au 25 octobre 2015

Pour aider (?), je mets un lien vers la température en France  dans quelques villes, paris, toulouse, grenoble, nantes, dunkerque (via ftp.ncdc.noaa.gov). Il sera possible d’utiliser la vraie température de 2015 pour faire les prévisions de consommation électrique.  Idéalement, je veux un fichier html (markdown) et 4 fichiers csv avec les 4 séries de prévision.

Moyennes de Prédicteurs

La dernière fois, nous avions vu 4 modèles de séries temporelles, pour modéliser la série hebdomadaire de températures. Le premier était obtenu avec auto.arima, les deux étaient ensuite ses modèles SARIMA (le second était une forte autocorélation à l’ordre 52, et le troisième avait une racine unité saisonnière), et le dernier était obtenu par la méthode de Buys-Ballot (cyclique, avec un bruit autocorrélé).

Au lieu de se demander lequel des quatre est le “meilleur”, on peut tenter de trouver une moyenne pondérée, \hat{y}_t = \sum_{i}\omega_i \hat{y}_t^{(i)}\omega=(\omega_i) est un point du simplexe, au sens où \omega_i\geq 0 et \omega_1+...+\omega_k=1. On veut ensuite trouver les “meilleurs” pondérations, par exemple au sens de la norme L^2. On va calibrer nos quatre modèles sur les m premières valeurs, et on compare ensuite la combinaison des prévions sur les 111 prochaines valeurs (et les vraies valeurs), \omega^\star=\text{argmin}\{\sum_{h=1}^{111}\left(y_{t+h}- \sum_{i}\omega_i \hat{y}_{t+h}^{(i)}\right)^2\}

Par exemple, on peut utiliser les 200 premières valeurs.

m=200
Ym=ts(y[1:m],start=1996,frequency=52)

On va ensuite estimer les 4 modèles sur ces 200 valeurs

model1=arima(Ym,order = c(4,0,0), 
seasonal = list(order = c(1,0,0))) 
model2=arima(Ym,order = c(1,0,0), 
seasonal = list(order = c(1,1,0)))
model3=arima(Ym,order = c(0,0,2), 
seasonal = list(order = c(0,1,1)))
regL=lm(Temp~0+as.factor(NumWeek),
data=electricite[1:m,])
E=residuals(regL)
model4=arima(ts(E,start=1996,
frequency=52),order = c(1,0,0))

On peut ensuite collecter les prédictions

DONNEES = data.frame(
  y1=predict(model1,n.ahead = 111)$pred,
  y2=predict(model2,n.ahead = 111)$pred,
  y3=predict(model3,n.ahead = 111)$pred,
  y4=predict(model4,n.ahead = 111)$pred+
     predict(regL,newdata=electricite[m+(1:111),]),
  vraie=electricite[m+(1:111),"Temp"])

Pour créer une combinaison linéaire des prévisions, on utilise

a=rep(1/4,4)
y_prev = as.matrix(DONNEES[,1:4]) %*% a

On peut ainsi visualiser les 4 prévisions, leur combinaison linéaire (ici avec des équi-pondérations), et leurs observations

plot(1:111,DONNEES[,5],type="l",lwd=2)
for(i in 1:4) lines(1:111,DONNEES[,i],col="blue")
lines(1:111,yp,col="red")

Pour trouver la valeur “optimale” des pondérations, qui minimisera la somme des carrés des erreurs, on utilise le code suivant (on fera une transformation des données, ici sur le simplexe, pour s’assurer que l’optimisation donnera un élément du simplexe)

erreur=function(a) sum((DONNEES[,1:4] 
%*% a-DONNEES[,5])^2)
erreur2=function(b){
  a=rep(NA,length(b)+1)
  a[1:length(b)]=exp(b)/(1+sum(exp(b)))
  a[length(b)+1]=1/(1+sum(exp(b)))
  erreur(a)
}

On obtient l’optimum, et on peut obtenir une prévision

b=optim(par=c(0,0,0),erreur2)$par
a=rep(NA,length(b)+1)
a[1:length(b)]=exp(b)/(1+sum(exp(b)))
a[length(b)+1]=1/(1+sum(exp(b)))
yp=DONNEES[,1:4] %*% a
lines(1:111,yp,col="red",lwd=4)

On peut alors faire changer le nombre d’observations utilisées

Techniquement, ca pose un petit soucis, car il faut s’assurer la convergence de deux algorithmes : celui d’estimation des paramètres des SARIMA, celui de recherche du paramètre de poids.

poids_optimaux=function(m){
  Ym=ts(y[1:m],start=1996,frequency=52)
  TRY=try(model1 <- arima(Ym,order = c(4,0,0), 
    seasonal = list(order = c(1,0,0))),FALSE)
  if(inherits(TRY, "try-error")) model1 <- arima(Ym,order = c(4,0,0), 
    seasonal = list(order = c(1,0,0)),method="CSS")
  TRY=try(model2 <- arima(Ym,order = c(1,0,0), 
    seasonal = list(order = c(1,1,0))),FALSE)
  if(inherits(TRY, "try-error")) model2 <- arima(Ym,order = c(1,0,0), 
    seasonal = list(order = c(1,1,0)),method="CSS")
  TRY=try(model3 <- arima(Ym,order = c(0,0,2), 
    seasonal = list(order = c(0,1,1))),FALSE)
  if(inherits(TRY, "try-error")) model3 <- arima(Ym,order = c(0,0,2), 
    seasonal = list(order = c(0,1,1)),method="CSS")
regL=lm(Temp~0+as.factor(NumWeek),data=electricite[1:m,])
  E=residuals(regL)
model4=arima(ts(E,start=1996,frequency=52),order = c(1,0,0))
  DONNEES = data.frame(
    y1=predict(model1,n.ahead = 111)$pred,
    y2=predict(model2,n.ahead = 111)$pred,
    y3=predict(model3,n.ahead = 111)$pred,
    y4=predict(model4,n.ahead = 111)$pred+
      predict(regL,newdata=electricite[m+(1:111),]),
    vraie=electricite[m+(1:111),"Temp"])
  erreur=function(a) sum((as.matrix(DONNEES[,1:4]) %*% a-DONNEES[,5])^2)
  erreur2=function(b){
    a=rep(NA,length(b)+1)
    a[1:length(b)]=exp(b)/(1+sum(exp(b)))
    a[length(b)+1]=1/(1+sum(exp(b)))
    erreur(a)
  }
  TRY=try(b <- optim(par=c(0,0,0),erreur2)$par, FALSE)
  a=rep(NA,length(b)+1)
  if(!inherits(TRY, "try-error")) {
  a[1:length(b)]=exp(b)/(1+sum(exp(b)))
  a[length(b)+1]=1/(1+sum(exp(b)))
  }
  return(a)
}

On peut alors regarder l’évolution des pondérations au cours du temps. On peut vectoriser, mais si l’algorithme aboutit sur une erreur, c’est plus simple d’identifier

VM=seq(200,480,by=10)
VA=matrix(NA,length(VM),4)
for(i in 1:length(VM)) VA[i,]=poids_optimaux(VM[i])
colnames(VA)=c("auto","SARMA","SARIMA","Buys Ballot")

(pour débugger, c’est plus simple les boucles). On peut alors obtenir le graphique suivant, avec en rose Buys-Ballot, et en vert le modèle SARIMA, avec une racine unité saisonnière,

tva=t(VA)
barplot(tva,legend = rownames(counts),
 col=c("yellow","orange","green","pink"))

On se rend compte qu’assez vite, le modèle le plus intéressant est le modèle de Buys Ballot.

Il est possible de changer de fonction de perte, par exemple si on s’intéresse à un quantile à 90%,

tau=.9
loss=function(e) e*(tau-(e<=0)*1)

Dans la fonction, on utilise alors

erreur=function(a) sum(loss(as.matrix(
DONNEES[,1:4]) %*% a-DONNEES[,5]))

C’est, cette fois, beaucoup plus incertain, même si les deux modèles les plus intéressants sont toujours le SARIMA et Buys-Ballot.

Consommation et production électrique

Une version temporaire des slides est en ligne. Pour le TP, on continuera à utiliser les données suivante

url="http://freakonometrics.free.fr/elec.RData"
download.file(url,"elec.RData")
load("elec.RData")
k=620
n=nrow(electricite)
passe=1:k
futur=(k+1):n
y=electricite$Load[passe]
plot(y,type="l")

La semaine passée, nous avions commencé de modéliser la série de température (cette dernière ayant un impact important comme déterminant de la consommation électrique)

y=electricite$Temp[passe]
plot(y,type="l")
Y=ts(y,start=1996,frequency=52)
library(forecast)

plot(y[(53:k)-52],y[53:k])
abline(lm(y[53:k]~y[(53:k)-52]),col="red")
abline(a=0,b=1,lty=2)

La série temporelle est relativement autocorrélée, à l’ordre 52

acf(y,lag=120)

model1=auto.arima(Y)
acf(residuals(model1),120)

On va garder ce modèle en mémoire, et regarder par la suite ses prévisions. Tentons ici un SARIMA, à la main

model2=arima(Y,order = c(0,0,0), 
seasonal = list(order = c(1,0,0)))
summary(model2)
model2=arima(Y,order = c(0,0,0), 
seasonal = list(order = c(2,0,0)))
summary(model2)
acf(residuals(model2),120)
model2=arima(Y,order = c(2,0,0), 
seasonal = list(order = c(1,0,0)))            
model2=arima(Y,order = c(1,0,0), 
seasonal = list(order = c(1,1,0)))
summary(model2)
acf(residuals(model2),120)

Tentons ensuite avec une racine unité saisonnière

Z=diff(Y,52)
model3=arima(Z,order = c(0,0,1), 
seasonal = list(order = c(0,0,1)))
summary(model3)
model3=arima(Z,order = c(0,0,2), 
seasonal = list(order = c(0,0,1)))
summary(model3)
acf(residuals(model3),120)
model3=arima(Y,order = c(0,0,2), 
seasonal = list(order = c(0,1,1)))
summary(model3)

On peut ensuite tenter un modèle à la Buys-Ballot

reg=lm(Temp~0+as.factor(NumWeek),
data=electricite[passe,])
cycle=predict(reg,newdata=electricite[passe,])
plot(y,type="l")
lines(cycle,col="red")
E=residuals(reg)
plot(E,type="l")
acf(E,120)
model4=arima(E,order = c(1,0,0))
summary(model4)

horizon=n-k
prev4E=predict(model4,n.ahead = horizon)$pred
prev4C=predict(reg,newdata=electricite[futur,])
prev4 =prev4E+prev4C
plot(y,type="l",xlim=c(0,n))
abline(v=k,col="red")
lines(prev4,col="blue")

prev1=predict(model1,n.ahead = horizon)$pred
plot(y,type="l",xlim=c(0,n))
abline(v=k,col="red")
lines(futur,prev1,col="purple")

prev2=predict(model2,n.ahead = horizon)$pred
plot(y,type="l",xlim=c(0,n))
abline(v=k,col="red")
lines(futur,prev2,col="green")

prev3=predict(model3,n.ahead = horizon)$pred
plot(y,type="l",xlim=c(0,n))
abline(v=k,col="red")
lines(futur,prev3,col="orange")

MODELES = data.frame(y1=prev1,
y2=prev2,y3=prev3,y4=prev4, 
y =electricite$Temp[futur])
plot(MODELES$y,type="l")
lines(1:111,MODELES$y1,col="purple")
lines(1:111,MODELES$y2,col="green")
lines(1:111,MODELES$y3,col="orange")
lines(1:111,MODELES$y4,col="blue")