Multi-Step Time Series Prediction with LSTM

Household Power Consumption Prediction using RNN-LSTM

A power outage can cause a huge economic loss. So it is very important to predict power use.

Smart meters and solar panels are now common. This means we have a lot of electricity usage data to work with.

Problem Statement :

Given that power consumption data for the previous week, we have to predict the power consumption for the next week.

Watch Full Video: 

Download dataset:

Download household_power_consumption.zip

Details:

UCI Household Electric Power Consumption dataset

Dataset Description:

The data was collected from December 2006 to November 2010. Power use in the household was recorded every minute.

It is a multivariate series comprised of seven variables

• global_active_power: The total active power consumed by the household (kilowatts).
• global_reactive_power: The total reactive power consumed by the household (kilowatts).
• voltage: Average voltage (volts).
• global_intensity: Average current intensity (amps).
• sub_metering_1: Active energy for kitchen (watt-hours of active energy).
• sub_metering_2: Active energy for laundry (watt-hours of active energy).
• sub_metering_3: Active energy for climate control systems (watt-hours of active energy).

This data is a multivariate time series of power values. We can use it to model and forecast future electricity use.

In this blog, we will build an encoder-decoder LSTM in TensorFlow to forecast household power consumption seven days ahead using the UCI household electric power dataset. LSTM networks handle long-range time dependencies through gated memory cells, which makes them a good fit for multi-step time-series forecasting.

Importing Libraries

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from numpy import nan

from tensorflow.keras import Sequential
from tensorflow.keras.layers import LSTM, Dense

from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import MinMaxScaler
python
#Reading the dataset

data = pd.read_csv('household_power_consumption.txt', sep = ';',
                  parse_dates = True,
                  low_memory = False)

#printing top rows

data.head()

#concatenating the date and time columns to 'date_time' columns

data['date_time'] = data['Date'].str.cat(data['Time'], sep= ' ')
data.drop(['Date', 'Time'], inplace= True, axis = 1)
data.head()

data.set_index(['date_time'], inplace=True)
data.head()

Next, we can mark all missing values indicated with a '?' character with a NaN value, which is a float.

#replacing each '?'characters with NaN value

data.replace('?', nan, inplace=True)

# cast all columns to float64 for uniform numeric computation

data = data.astype('float')

#information of the dataset

data.info()

Index: 2075259 entries, 16/12/2006 17:24:00 to 26/11/2010 21:02:00
Data columns (total 7 columns):
Global_active_power      float64
Global_reactive_power    float64
Voltage                  float64
Global_intensity         float64
Sub_metering_1           float64
Sub_metering_2           float64
Sub_metering_3           float64
dtypes: float64(7)
memory usage: 126.7+ MB

#checking the null values

np.isnan(data).sum()

Global_active_power      25979
Global_reactive_power    25979
Voltage                  25979
Global_intensity         25979
Sub_metering_1           25979
Sub_metering_2           25979
Sub_metering_3           25979
dtype: int64

We also need to fill in the missing values now that they have been marked.

A very simple approach would be to copy the observation from the same time the day before. We can implement this in a function named fill_missing() that will take the NumPy array of the data and copy values from exactly 24 hours ago.

def fill_missing(data):
    one_day = 24*60
    for row in range(data.shape[0]):
        for col in range(data.shape[1]):
            if np.isnan(data[row, col]):
                data[row, col] = data[row-one_day, col]

fill_missing(data.values)

#checking the nan values

np.isnan(data).sum()

Global_active_power      0
Global_reactive_power    0
Voltage                  0
Global_intensity         0
Sub_metering_1           0
Sub_metering_2           0
Sub_metering_3           0
dtype: int64

data.info()

Index: 2075259 entries, 16/12/2006 17:24:00 to 26/11/2010 21:02:00
Data columns (total 7 columns):
Global_active_power      float64
Global_reactive_power    float64
Voltage                  float64
Global_intensity         float64
Sub_metering_1           float64
Sub_metering_2           float64
Sub_metering_3           float64
dtypes: float64(7)
memory usage: 126.7+ MB

#printing the shape of the data

data.shape

(2075259, 7)

Here, we can observe that we have 2075259 datapoints and 7 features

data.head()

Prepare power consumption for each day

We can now save the cleaned-up version of the dataset to a new file; in this case we will just change the file extension to .csv and save the dataset as 'cleaned_data.csv'.

#conversion of dataframe to .csv

data.to_csv('cleaned_data.csv')

#reading the dataset

dataset = pd.read_csv('cleaned_data.csv', parse_dates = True, index_col = 'date_time', low_memory = False)

#printing the top rows

dataset.head()

#printing the bottom rows

dataset.tail()

Exploratory Data Analysis

# resample to daily bins, summing all timestamps within each day

data = dataset.resample('D').sum()

#data after sampling it into daywise manner

data.head()

Plotting the all features in various time stamps

fig, ax = plt.subplots(figsize=(18,18))

for i in range(len(data.columns)):
    plt.subplot(len(data.columns), 1, i+1)
    name = data.columns[i]
    plt.plot(data[name])
    plt.title(name, y=0, loc = 'right')
    plt.yticks([])
plt.show()
fig.tight_layout()

Exploring Active power consumption for each year

# four full years of data available

years = ['2007', '2008', '2009', '2010']

Year wise plotting of feature Global_active_power

fig, ax = plt.subplots(figsize=(18,18))

for i in range(len(years)):
    plt.subplot(len(years), 1, i+1)
    year = years[i]
    active_power_data = data[str(year)]
    plt.plot(active_power_data['Global_active_power'])
    plt.title(str(year), y = 0, loc = 'left')
plt.show()
fig.tight_layout()

#for year 2006

data['2006']

Power consumption distribution with histogram

Year wise histogram plot of feature Global_active_power

fig, ax = plt.subplots(figsize=(18,18))

for i in range(len(years)):
    plt.subplot(len(years), 1, i+1)
    year = years[i]
    active_power_data = data[str(year)]
    active_power_data['Global_active_power'].hist(bins = 200)
    plt.title(str(year), y = 0, loc = 'left')
plt.show()
fig.tight_layout()

Histogram plot for All Features

fig, ax = plt.subplots(figsize=(18,18))

for i in range(len(data.columns)):
    plt.subplot(len(data.columns), 1, i+1)
    name = data.columns[i]
    data[name].hist(bins=200)
    plt.title(name, y=0, loc = 'right')
    plt.yticks([])
plt.show()
fig.tight_layout()

Plot power consumption hist for each month of 2007

months = [i for i in range(1,13)]

fig, ax = plt.subplots(figsize=(18,18))

for i in range(len(months)):
    ax = plt.subplot(len(months), 1, i+1)
    month = '2007-' + str(months[i])
    active_power_data = dataset[month]
    active_power_data['Global_active_power'].hist(bins = 100)
    ax.set_xlim(0,5)
    plt.title(month, y = 0, loc = 'right')
plt.show()
fig.tight_layout()

Observation :

1. From the above diagram we can say that power consumption in the month of Nov, Dec, Jan, Feb, Mar is more as there is a long tail as compare to other months.

2. It also shows that the during the winter seasons, the heating systems are used and not in summer.

3. The above graph is highly concentrated on 0.3W and 1.3W.

Active Power Uses Prediction

What can we predict

• Forecast hourly consumption for the next day.
• Forecast daily consumption for the next week.
• Forecast daily consumption for the next month.
• Forecast monthly consumption for the next year.

Modeling Methods

There are many modeling methods and few of those are as follows

• Naive Methods -> Naive methods would include methods that make very simple, but often very effective assumptions.
• Classical Linear Methods -> Classical linear methods include techniques are very effective for univariate time series forecasting
• Machine Learning Methods -> Machine learning methods require that the problem be framed as a supervised learning problem.K-nearest neighbors.
• SVM
• Decision trees
• Random forest
• Gradient boosting machines

#top rows

data.head()

#printing last rows

data.tail()

# train: up to end of 2009; test: 2010 onwards

data_train = data.loc[:'2009-12-31', :]['Global_active_power']
data_train.head()

date_time
2006-12-16    1209.176
2006-12-17    3390.460
2006-12-18    2203.826
2006-12-19    1666.194
2006-12-20    2225.748
Freq: D, Name: Global_active_power, dtype: float64

data_test = data['2010']['Global_active_power']
data_test.head()

date_time
2010-01-01    1224.252
2010-01-02    1693.778
2010-01-03    1298.728
2010-01-04    1687.440
2010-01-05    1320.158
Freq: D, Name: Global_active_power, dtype: float64

data_train.shape

(1112,)

data_test.shape

(345,)

#training data

data_train.head(14)

date_time
2006-12-16    1209.176
2006-12-17    3390.460
2006-12-18    2203.826
2006-12-19    1666.194
2006-12-20    2225.748
2006-12-21    1723.288
2006-12-22    2341.338
2006-12-23    4773.386
2006-12-24    2550.012
2006-12-25    2743.120
2006-12-26    3934.110
2006-12-27    1528.760
2006-12-28    2072.638
2006-12-29    3174.392
Freq: D, Name: Global_active_power, dtype: float64

#converting the data into numpy array

data_train = np.array(data_train)

# split data into weekly windows of 7 days

X_train, y_train = [], []

for i in range(7, len(data_train)-7):
    X_train.append(data_train[i-7:i])
    y_train.append(data_train[i:i+7])

#converting list to numpy array

X_train, y_train = np.array(X_train), np.array(y_train)

#shape of train and test dataset

X_train.shape, y_train.shape

((1098, 7), (1098, 7))

#printing the ytrain value

pd.DataFrame(y_train).head()

#Normalising the dataset between 0 and 1

x_scaler = MinMaxScaler()
X_train = x_scaler.fit_transform(X_train)

#Normalising the dataset

y_scaler = MinMaxScaler()
y_train = y_scaler.fit_transform(y_train)

pd.DataFrame(X_train).head()

#converting to 3 dimension

X_train = X_train.reshape(1098, 7, 1)

X_train.shape

(1098, 7, 1)

#building sequential model using Keras

reg = Sequential()
reg.add(LSTM(units = 200, activation = 'relu', input_shape=(7,1)))
reg.add(Dense(7))

# MSE loss, Adam optimizer

reg.compile(loss='mse', optimizer='adam')

#training the model

reg.fit(X_train, y_train, epochs = 100)

Train on 1098 samples
Epoch 1/100
1098/1098 [==============================] - 2s 2ms/sample - loss: 0.0626
Epoch 2/100
1098/1098 [==============================] - 0s 296us/sample -
.
.
.
.
.
Epoch 99/100
1098/1098 [==============================] - 0s 270us/sample - loss: 0.0228
Epoch 100/100
1098/1098 [==============================] - 0s 269us/sample - loss: 0.0228

#testing dataset

data_test = np.array(data_test)

# split test data into weekly windows of 7 days

X_test, y_test = [], []

for i in range(7, len(data_test)-7):
    X_test.append(data_test[i-7:i])
    y_test.append(data_test[i:i+7])

X_test, y_test = np.array(X_test), np.array(y_test)

X_test = x_scaler.transform(X_test)
y_test = y_scaler.transform(y_test)

#converting to 3 dimension

X_test = X_test.reshape(331,7,1)

X_test.shape

(331, 7, 1)

y_pred = reg.predict(X_test)

#bringing y_pred values to their original forms by using inverse transform

y_pred = y_scaler.inverse_transform(y_pred)

y_pred

array([[1508.9413 , 1476.1537 , 1487.5676 , ..., 1484.8464 , 1459.3864 ,
        1551.5675 ],
       [1158.2788 , 1287.0326 , 1346.428  , ..., 1430.5685 , 1420.6346 ,
        1472.5759 ],
       [1571.7665 , 1507.0337 , 1516.5574 , ..., 1432.5813 , 1393.9161 ,
        1504.1714 ],
       ...,
       [ 952.85785,  852.4236 ,  933.62585, ...,  800.12006,  831.2844 ,
        1005.20844],
       [1579.4896 , 1353.6078 , 1278.9501 , ...,  981.4198 ,  967.6466 ,
        1146.7898 ],
       [1629.0509 , 1392.7751 , 1288.7218 , ..., 1052.977  , 1070.8586 ,
        1243.1346 ]], dtype=float32)

y_true = y_scaler.inverse_transform(y_test)

y_true

array([[ 555.664, 1593.318, 1504.82 , ...,    0.   , 1995.796, 2116.224],
       [1593.318, 1504.82 , 1383.18 , ..., 1995.796, 2116.224, 2196.76 ],
       [1504.82 , 1383.18 ,    0.   , ..., 2116.224, 2196.76 , 2150.112],
       ...,
       [1892.998, 1645.424, 1439.426, ..., 1973.382, 1109.574,  529.698],
       [1645.424, 1439.426, 2035.418, ..., 1109.574,  529.698, 1612.092],
       [1439.426, 2035.418, 1973.382, ...,  529.698, 1612.092, 1579.692]])

def evaluate_model(y_true, y_predicted):
    scores = []

    #calculate scores for each day
    for i in range(y_true.shape[1]):
        mse = mean_squared_error(y_true[:, i], y_predicted[:, i])
        rmse = np.sqrt(mse)
        scores.append(rmse)

    #calculate score for whole prediction
    total_score = 0
    for row in range(y_true.shape[0]):
        for col in range(y_predicted.shape[1]):
            total_score = total_score + (y_true[row, col] - y_predicted[row, col])**2
    total_score = np.sqrt(total_score/(y_true.shape[0]*y_predicted.shape[1]))

    return total_score, scores

evaluate_model(y_true, y_pred)

(579.2827596682928, [598.0411885086157, 592.5770673397814, 576.1153945912635, 563.9396525162248, 576.5479538079353, 570.7699415990154, 576.2430188855649])

#standard deviation

np.std(y_true[0])

710.0253857243853