Multivariable Time Series Forecasting with Neural Networks

Hello, today we are going to create a model to predict a forecast from a dataset from Intellisense.io.

IntelliSense.io a Market Leader in the Industrial Internet of Things and Artificial Intelligence .IntelliSense.io a Cambridge (UK) based provider of Optimisation as a Service (OaaS) applications for the Natural Resources Industries.

The company’s applications are delivered from their Industrial Internet of Things (IIoT) and Artificial Intelligence (AI) powered real-time decision-making platform (Brains.app), through a No Capital Expenditure model.

These applications have been deployed in the Mining industry across Latin America and Kazakhstan for base and precious metal mining operations. IntelliSense.io is a privately-held company based in Cambridge, UK, with offices in Almaty, Kazakhstan, Barcelona, Spain and Santiago, Chile.

Introduction

The goal of the project is to devise and implement a prototype data-driven prediction model of a Semi-autogenous (SAG) Mill.

Semi-autogenous (SAG) mills are central pieces of equipment for large mining operations. They are commonly used in the secondary crushing stage to break down larger rocks from the pit for further processing. As a very rough analogy, a SAG mill works like a giant washing machine for rocks and steel balls.

Description of the problem

IntelliSense.io provides one year of data based on real measurements from a client SAG mill.

Types of data tracked include performance variables, which are used to monitor the operation of the mill.

We are interested to predict:

1. Power Draw (MW) — Power drawn by the mill motor
2. Bearing Pressure (kPa) — Pressure on the mill supports. Can be thought of as the weight of the mill.

By using the data as an input:

1. Speed (RPM) — Rotation speed of the mill
2. Conveyor Belt Feed Rate (t/h) — Mass flow rate of material into the mill
3. Dilution Flow Rate (m³/h) — Volume flow rate of water into the mill
4. Conveyor Belt Fines (%) — Percent of material in the feed classified as “fine”

For this project I need the libraries:

from numpy import concatenate
from matplotlib import pyplot
from pandas import read_csv
from pandas import DataFrame
from pandas import concat
from sklearn.preprocessing import MinMaxScaler
from sklearn.preprocessing import LabelEncoder
from sklearn.metrics import mean_squared_error
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import LSTM

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

import scipy.stats as stats

from pandas import read_csv
from datetime import datetime
import pandas as pd
from pandas import read_csv
from matplotlib import pyplot
from pandas import read_csv
from matplotlib import pyplot

Load and inspect data

data= pd.read_csv("data/sag_data_train.csv", index_col="Time", parse_dates=True)

data.head()

	Bearing Pressure (kPa)	Power Draw (MW)	Speed (RPM)	Dilution Flow Rate (m3/h)	Conveyor Belt Feed Rate (t/h)	Conveyor Belt PSD Fines (%)
Time
2015-09-15 00:00:00	5488.175540	11.737357	7.843532	1030.590108	2935.660276	38.641018
2015-09-15 00:01:00	5423.930126	11.543755	7.681607	1039.869847	2928.333772	45.243656
2015-09-15 00:02:00	5502.058523	11.169525	7.514173	1033.237205	2919.128115	38.716221
2015-09-15 00:03:00	5477.862749	11.035091	7.592248	1035.075573	2985.500811	42.860703
2015-09-15 00:04:00	5508.013201	11.418827	7.784895	1042.189406	2905.052105	50.524544

Before develop any model of machine learning first we analyze the data Lets find if there is any liner correlation between variables

corr = data.corr()
fig = plt.figure()
ax = fig.add_subplot(111)
cax = ax.matshow(corr,cmap='coolwarm', vmin=-1, vmax=1)
fig.colorbar(cax)
ticks = np.arange(0,len(data.columns),1)
ax.set_xticks(ticks)
plt.xticks(rotation=90)
ax.set_yticks(ticks)
ax.set_xticklabels(data.columns)
ax.set_yticklabels(data.columns)
plt.show()

def personserie(df,paraX,paraY):
    df = df[[paraX, paraY]]
    overall_pearson_r = df.corr().iloc[0,1]
    print(f"Pandas computed Pearson r: {overall_pearson_r}")
    r, p = stats.pearsonr(df.dropna()[paraX], df.dropna()[paraY])
    print(f"Scipy computed Pearson r: {r} and p-value: {p}")
    f,ax=plt.subplots(figsize=(14,3))
    df.rolling(window=30,center=True).median().plot(ax=ax)
    
    ax.set(xlabel='Frame',ylabel='Smiling evidence',title=f"Overall Pearson r = {np.round(overall_pearson_r,2)}");
       

personserie(data,'Power Draw (MW)','Speed (RPM)')

Pandas computed Pearson r: 0.9804985064390073
Scipy computed Pearson r: 0.9804985064390568 and p-value: 0.0

This means that best option to predict the PowerDraw should be use a linear regresion and Bearing Pressure by times series forecasting method. The project will be splitted into the following sections:

Section 1 - Bearing Pressure

Section 2 - PowerDraw

Section 1 - Bearing Pressure

Time series forecasting problems must be re-framed as supervised learning problems. From a sequence to pairs of input and output sequences.Technically, in time series forecasting terminology the current time (t) and future times (t+1, t+n) are forecast times and past observations (t-1, t-n) are used to make forecasts. Positive and negative shifts can be used to create a new DataFrame from a time series with sequences of input and output patterns for a supervised learning problem. Further, the shift function also works on so-called multivariate time series problems. All variates in the time series can be shifted forward or backward to create multivariate input and output sequences.

We will define a new Python function named series_to_supervised() that takes a univariate or multivariate time series and frames it as a supervised learning dataset. The function takes four arguments:

data: Sequence of observations as a list or 2D NumPy array. Required.

n_in: Number of lag observations as input (X). Values may be between [1..len(data)] Optional. Defaults to 1.

n_out: Number of observations as output (y). Values may be between [0..len(data)-1]. Optional. Defaults to 1.

dropnan: Boolean whether or not to drop rows with NaN values. Optional. Defaults to True. The function returns a single value:

return: Pandas DataFrame of series framed for supervised learning.

# convert series to supervised learning
def series_to_supervised(data, n_in=1, n_out=1, dropnan=True):
    n_vars = 1 if type(data) is list else data.shape[1]
    df = DataFrame(data)
    cols, names = list(), list()
    # input sequence (t-n, ... t-1)
    for i in range(n_in, 0, -1):
        cols.append(df.shift(i))
        names += [('var%d(t-%d)' % (j+1, i)) for j in range(n_vars)]
    # forecast sequence (t, t+1, ... t+n)
    for i in range(0, n_out):
        cols.append(df.shift(-i))
        if i == 0:
            names += [('var%d(t)' % (j+1)) for j in range(n_vars)]
        else:
            names += [('var%d(t+%d)' % (j+1, i)) for j in range(n_vars)]
    # put it all together
    agg = concat(cols, axis=1)
    agg.columns = names
    # drop rows with NaN values
    if dropnan:
        agg.dropna(inplace=True)
    return agg

# load dataset
df_train = pd.read_csv("data/sag_data_train.csv", index_col="Time", parse_dates=True)

#df_train.head()

DATA PREPARATION

The first step is to prepare the Bearing Pressure (kPa) dataset. This involves framing the dataset as a supervised learning problem and normalizing the input variables. We will frame the supervised learning problem as predicting the Bearing Pressure at the current time (t) given the Speed (RPM) Dilution Flow Rate (m3/h) Conveyor Belt Feed Rate (t/h) Conveyor Belt PSD Fines (%) measurements at the prior time step.

dataset=df_train
values = df_train.values

# specify columns to plot
groups = [0, 1, 2, 3, 4, 5]
i = 1
# plot each column
pyplot.figure()
for group in groups:
    pyplot.subplot(len(groups), 1, i)
    pyplot.plot(values[:, group])
    pyplot.title(dataset.columns[group], y=0.5, loc='right')
    i += 1
pyplot.show()

The first step is to prepare the Bearing Pressure (kPa) dataset for the LSTM. This involves framing the dataset as a supervised learning problem and normalizing the input variables.We will frame the supervised learning problem as predicting the Bearing Pressure (kPa) at the current time (t) given the Power Draw (MW),Speed (RPM) Dilution Flow Rate (m3/h) Conveyor Belt Feed Rate (t/h) Conveyor Belt PSD Fines (%) parameters at the prior time step.

Fortunatelly we dont have Categorical Data on the features, so we dont need Convert Categorical Data to Numerical Data by Integer Encoding or One-Hot Encoding.

# integer encode direction
#encoder = LabelEncoder()
#values[:,number_col] = encoder.fit_transform(values[:,number_col])

# ensure all data is float
values = values.astype('float32')

Next, all features are normalized, then the dataset is transformed into a supervised learning problem.

# normalize features
scaler = MinMaxScaler(feature_range=(0, 1))
scaled = scaler.fit_transform(values)

scaled[0]

array([0.91485214, 0.7853231 , 0.8984538 , 0.62976223, 0.81679434,
       0.45599252], dtype=float32)

Multivariate Forecasting

This is where we may have observations of multiple different measures and an interest in forecasting one or more of them. For example, we may have two sets of time series observations obs1 and obs2 and we wish to forecast one or both of these. We can call series_to_supervised() in exactly the same way.

# frame as supervised learning
reframed = series_to_supervised(scaled, 1, 1)

printing the new framing of the data, showing an input pattern with one time step for both variables and an output pattern of one time step for both variables. Again, depending on the specifics of the problem, the division of columns into X and Y components can be chosen arbitrarily, such as if the current observation of var1 was also provided as input and only var2 was to be predicted.

#len(reframed)

reframed

	var1(t-1)	var2(t-1)	var3(t-1)	var4(t-1)	var5(t-1)	var6(t-1)	var1(t)	var2(t)	var3(t)	var4(t)	var5(t)	var6(t)
1	0.914852	0.785323	0.898454	0.629762	0.816794	0.455993	0.904073	0.772450	0.879906	0.635433	0.814756	0.533909
2	0.904073	0.772450	0.879906	0.635433	0.814756	0.533909	0.917181	0.747567	0.860727	0.631380	0.812195	0.456880
3	0.917181	0.747567	0.860727	0.631380	0.812195	0.456880	0.913122	0.738628	0.869670	0.632503	0.830662	0.505788
4	0.913122	0.738628	0.869670	0.632503	0.830662	0.505788	0.918181	0.764143	0.891737	0.636850	0.808278	0.596227
5	0.918181	0.764143	0.891737	0.636850	0.808278	0.596227	0.914700	0.782284	0.909612	0.631457	0.819966	0.558852
...	...	...	...	...	...	...	...	...	...	...	...	...
138236	0.888515	0.902305	0.994407	0.586262	0.794455	0.500200	0.881775	0.913288	0.997608	0.576388	0.800384	0.593150
138237	0.881775	0.913288	0.997608	0.576388	0.800384	0.593150	0.884448	0.909902	0.985473	0.596625	0.782906	0.530104
138238	0.884448	0.909902	0.985473	0.596625	0.782906	0.530104	0.876622	0.895403	0.962863	0.570117	0.800633	0.477679
138239	0.876622	0.895403	0.962863	0.570117	0.800633	0.477679	0.873892	0.876178	0.943823	0.581300	0.801972	0.565912
138240	0.873892	0.876178	0.943823	0.581300	0.801972	0.565912	0.880246	0.870881	0.934119	0.593477	0.796080	0.565529

138234 rows × 12 columns

The variables for the time to be predicted (t) are then removed.

It is printed the first 5 rows of the transformed dataset.

# drop columns we don't want to predict
reframed.drop(reframed.columns[[7,8,9,10,11]], axis=1, inplace=True)
print(reframed.head())

   var1(t-1)  var2(t-1)  var3(t-1)  var4(t-1)  var5(t-1)  var6(t-1)   var1(t)
1   0.914852   0.785323   0.898454   0.629762   0.816794   0.455993  0.904073
2   0.904073   0.772450   0.879906   0.635433   0.814756   0.533909  0.917181
3   0.917181   0.747567   0.860727   0.631380   0.812195   0.456880  0.913122
4   0.913122   0.738628   0.869670   0.632503   0.830662   0.505788  0.918181
5   0.918181   0.764143   0.891737   0.636850   0.808278   0.596227  0.914700

We can see the 6 input variables (input series) and the 1 output variable (Bearing Pressure (kPa) level at the current time).The Power Draw (MW) it is not needed to include as an input because it linear correlated with the Speed (RPM). However in order to develop the model we include them. Also can be removed the outliers. Additionally shold be possible make all series stationary with differencing and seasonal adjustment. But for the current times series we keep as it is. The scope for now is not optimization. Notice that we can choose the inputs in this part and try to identify which inputs are more important by using PCA but the current scope for now is show the method to perform this kind of predictions.

FITTING PROCEDURE

We will fit an LSTM on the multivariate input data. First, we must split the prepared dataset into train and test sets. we will only fit the model taking the first dataset sag_data_train.csv, the dataset ag_data_test.csv will be used only to test our future models. The standard procedure of splitting is given by the splits the dataset into train 80% of data and for the test the 20% of data, however for the reason we want predict 5 minutes in the future the testing will be considered 25 minutes by a window for the training test, We splits the train and test sets into input and output variables.Finally, the inputs (X) are reshaped into the 3D format expected by LSTMs

Let us know what is the index time of the first row of the traing set

df_train.index[0]

Timestamp('2015-09-15 00:00:00')

and the final

df_train.index[-1]

Timestamp('2015-12-20 00:00:00')

#df_train.index[138240]

So the training set is given between 2015-09-15 -> 2015-12-20

#df_train.shape

#df_train.tail(1).index

#df_train['2015-09-15 00:00:00':].head(1).index

Let us count the number of days

from datetime import date
f_date = date(2015,9,15)
l_date = date(2015,12,20)
delta = l_date - f_date
print('Between the initial date and final day there are',delta.days,'days')

Between the initial date and final day there are 96 days

print('Which converted to minutes is ',delta.days*1440)
print( 'which is ',len(df_train)-1,'rows in our dataframe train')

Which converted to minutes is  138240
which is  138240 rows in our dataframe test

The stantard procedure take the 80% of training and 20% of testing

print('the training will have',int(delta.days*80/100*1440),'minutes')
print('the test will have',int(delta.days*20/100*1440),'minutes')

the training will have 110592 minutes
the test will have 27648 minutes

# split into train and test sets
values = reframed.values
n_train_minutes = int(delta.days*80/100*1440)
train = values[:n_train_minutes, :]
#test = values[n_train_minutes:, :]

However I would like take a windows to test of 25 min

window=25
test = values[n_train_minutes:n_train_minutes+window, :]

len(test)

25

# split into input and outputs
train_X, train_y = train[:, :-1], train[:, -1]
test_X, test_y = test[:, :-1], test[:, -1]

# reshape input to be 3D [samples, timesteps, features]
train_X = train_X.reshape((train_X.shape[0], 1, train_X.shape[1]))
test_X = test_X.reshape((test_X.shape[0], 1, test_X.shape[1]))
print(train_X.shape, train_y.shape, test_X.shape, test_y.shape)

(110592, 1, 6) (110592,) (25, 1, 6) (25,)

print('The shape of the train and test input and output sets with about', train_X.shape[0],' minutes of data for training and about',  test_X.shape[0],' minutes for testing, which is great.')

The shape of the train and test input and output sets with about 110592  minutes of data for training and about 25  minutes for testing, which is great.

LSTM model

We will define the LSTM with 50 neurons in the first hidden layer and 1 neuron in the output layer for predicting Bearing Pressure. The input shape will be 1 time step with 8 features.We will use the Mean Absolute Error (MAE) loss function and the efficient Adam version of stochastic gradient descent.The model will be fit for 50 training epochs with a batch size of 72.

from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense, LSTM, Dropout, RepeatVector, TimeDistributed

# design network
model = Sequential()
model.add(LSTM(50, input_shape=(train_X.shape[1], train_X.shape[2])))
model.add(Dense(1))
model.compile(loss='mae', optimizer='adam')
model.summary()

Model: "sequential"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
lstm (LSTM)                  (None, 50)                11400     
_________________________________________________________________
dense (Dense)                (None, 1)                 51        
=================================================================
Total params: 11,451
Trainable params: 11,451
Non-trainable params: 0
_________________________________________________________________

#from keras.utils.vis_utils import plot_model
#plot_model(model, to_file='model_plot.png', show_shapes=True, show_layer_names=True)

Finally, we keep track of both the training and test loss during training by setting the validation_data argument in the fit() function.

# fit network
history = model.fit(train_X, train_y, epochs=50, batch_size=72, validation_data=(test_X, test_y), verbose=2, shuffle=False)

Epoch 1/50
1536/1536 - 2s - loss: 0.0479 - val_loss: 0.2239
Epoch 2/50
1536/1536 - 1s - loss: 0.0131 - val_loss: 0.0215
Epoch 3/50
1536/1536 - 1s - loss: 0.0110 - val_loss: 0.0024
Epoch 4/50
1536/1536 - 1s - loss: 0.0108 - val_loss: 8.6463e-04
Epoch 5/50
1536/1536 - 1s - loss: 0.0105 - val_loss: 0.0019
Epoch 6/50
1536/1536 - 1s - loss: 0.0103 - val_loss: 0.0037
Epoch 7/50
1536/1536 - 1s - loss: 0.0106 - val_loss: 0.0029
Epoch 8/50
1536/1536 - 1s - loss: 0.0103 - val_loss: 0.0023
Epoch 9/50
1536/1536 - 1s - loss: 0.0101 - val_loss: 7.4642e-04
Epoch 10/50
1536/1536 - 1s - loss: 0.0101 - val_loss: 5.3502e-04
Epoch 11/50
1536/1536 - 1s - loss: 0.0098 - val_loss: 0.0019
Epoch 12/50
1536/1536 - 1s - loss: 0.0098 - val_loss: 0.0089
Epoch 13/50
1536/1536 - 1s - loss: 0.0096 - val_loss: 0.0077
Epoch 14/50
1536/1536 - 1s - loss: 0.0095 - val_loss: 0.0084
Epoch 15/50
1536/1536 - 1s - loss: 0.0094 - val_loss: 0.0080
Epoch 16/50
1536/1536 - 1s - loss: 0.0093 - val_loss: 0.0128
Epoch 17/50
1536/1536 - 1s - loss: 0.0092 - val_loss: 0.0051
Epoch 18/50
1536/1536 - 1s - loss: 0.0091 - val_loss: 0.0118
Epoch 19/50
1536/1536 - 1s - loss: 0.0090 - val_loss: 0.0126
Epoch 20/50
1536/1536 - 1s - loss: 0.0089 - val_loss: 0.0125
Epoch 21/50
1536/1536 - 1s - loss: 0.0089 - val_loss: 0.0016
Epoch 22/50
1536/1536 - 1s - loss: 0.0089 - val_loss: 0.0018
Epoch 23/50
1536/1536 - 1s - loss: 0.0088 - val_loss: 8.1201e-04
Epoch 24/50
1536/1536 - 1s - loss: 0.0087 - val_loss: 0.0117
Epoch 25/50
1536/1536 - 1s - loss: 0.0088 - val_loss: 0.0091
Epoch 26/50
1536/1536 - 1s - loss: 0.0089 - val_loss: 0.0090
Epoch 27/50
1536/1536 - 1s - loss: 0.0087 - val_loss: 0.0111
Epoch 28/50
1536/1536 - 1s - loss: 0.0085 - val_loss: 2.9275e-04
Epoch 29/50
1536/1536 - 1s - loss: 0.0086 - val_loss: 0.0062
Epoch 30/50
1536/1536 - 1s - loss: 0.0085 - val_loss: 0.0055
Epoch 31/50
1536/1536 - 1s - loss: 0.0085 - val_loss: 0.0099
Epoch 32/50
1536/1536 - 1s - loss: 0.0088 - val_loss: 0.0051
Epoch 33/50
1536/1536 - 1s - loss: 0.0084 - val_loss: 0.0102
Epoch 34/50
1536/1536 - 1s - loss: 0.0086 - val_loss: 0.0047
Epoch 35/50
1536/1536 - 1s - loss: 0.0086 - val_loss: 0.0131
Epoch 36/50
1536/1536 - 1s - loss: 0.0089 - val_loss: 6.6564e-04
Epoch 37/50
1536/1536 - 1s - loss: 0.0086 - val_loss: 0.0144
Epoch 38/50
1536/1536 - 1s - loss: 0.0089 - val_loss: 0.0100
Epoch 39/50
1536/1536 - 1s - loss: 0.0088 - val_loss: 0.0015
Epoch 40/50
1536/1536 - 1s - loss: 0.0086 - val_loss: 0.0081
Epoch 41/50
1536/1536 - 1s - loss: 0.0085 - val_loss: 0.0115
Epoch 42/50
1536/1536 - 1s - loss: 0.0087 - val_loss: 0.0102
Epoch 43/50
1536/1536 - 1s - loss: 0.0087 - val_loss: 0.0129
Epoch 44/50
1536/1536 - 1s - loss: 0.0086 - val_loss: 0.0030
Epoch 45/50
1536/1536 - 1s - loss: 0.0083 - val_loss: 0.0132
Epoch 46/50
1536/1536 - 1s - loss: 0.0085 - val_loss: 0.0043
Epoch 47/50
1536/1536 - 1s - loss: 0.0082 - val_loss: 0.0127
Epoch 48/50
1536/1536 - 1s - loss: 0.0083 - val_loss: 0.0077
Epoch 49/50
1536/1536 - 1s - loss: 0.0084 - val_loss: 0.0132
Epoch 50/50
1536/1536 - 1s - loss: 0.0083 - val_loss: 0.0062

At the end of the run both the training and test loss are plotted.

# plot history
pyplot.plot(history.history['loss'], label='train - loss')
pyplot.plot(history.history['val_loss'], label='test- val_loss')
pyplot.legend()
pyplot.show()

In machine learning and deep learning there are basically three cases:

1) Underfitting This is the only case where loss > validation_loss, but only slightly, if loss is far higher than validation_loss, please post your code and data so that we can have a look at

2) Overfitting loss « validation_loss This means that your model is fitting very nicely the training data but not at all the validation data, in other words it’s not generalizing correctly to unseen data

3) Perfect fitting loss == validation_loss If both values end up to be roughly the same and also if the values are converging (plot the loss over time) then chances are very high that you are doing it right

# make a prediction
yhat = model.predict(test_X)
test_X = test_X.reshape((test_X.shape[0], test_X.shape[2]))

# invert scaling for forecast
inv_yhat = concatenate((yhat, test_X[:, 1:]), axis=1)

inv_yhat = scaler.inverse_transform(inv_yhat)

inv_yhat = inv_yhat[:,0]

# invert scaling for actual
test_y = test_y.reshape((len(test_y), 1))
inv_y = concatenate((test_y, test_X[:, 1:]), axis=1)
inv_y = scaler.inverse_transform(inv_y)
inv_y = inv_y[:,0]

from math import sqrt

# calculate RMSE
rmse = sqrt(mean_squared_error(inv_y, inv_yhat))
print('Test RMSE: %.3f' % rmse)

Test RMSE: 37.242

Prediction

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
plt.style.use("seaborn")

With the model created lets take the test data provided and we try to perform a prediction with our working model

df_test = pd.read_csv("data/sag_data_test.csv", index_col="Time", parse_dates=True)

We want to predict the first 5 minutes of the test input

print( 'Starting from ',df_test.index.min(),'  plus 5 minutes')

Starting from  2015-12-21 00:00:00   plus 5 minutes

df_test_5min=df_test['2015-12-21 00:00:00':'2015-12-21 00:05:00']

#df_test_5min.plot(subplots=True, layout=(-1, 1), lw=1, figsize=(12, 6))
#plt.tight_layout()

In order to use our model, we have to makeup our input as the same way we used in our training datset

values_test=df_test_5min.values
values_test = values_test.astype('float32')
# normalize features
scaler = MinMaxScaler(feature_range=(0, 1))
scaled_test = scaler.fit_transform(values_test)
# frame as supervised learning
reframed_test = series_to_supervised(scaled_test, 1, 1)
# drop columns we don't want to predict
reframed_test.drop(reframed_test.columns[[7,8,9,10,11]], axis=1, inplace=True)
#  test sets
values_test = reframed_test.values
window=5
# reshape input to be 3D [samples, timesteps, features]
test = values_test[:window, :]
#len(test)
test_X, test_y = test[:, :-1], test[:, -1]
test_X = test_X.reshape((test_X.shape[0], 1, test_X.shape[1]))
yhat = model.predict(test_X)

test_X = test_X.reshape((test_X.shape[0], test_X.shape[2]))
# invert scaling for forecast
inv_yhat = concatenate((yhat, test_X[:, 1:]), axis=1)
inv_yhat = scaler.inverse_transform(inv_yhat)
inv_yhat = inv_yhat[:,0]
# invert scaling for actual
test_y = test_y.reshape((len(test_y), 1))
inv_y = concatenate((test_y, test_X[:, 1:]), axis=1)
inv_y = scaler.inverse_transform(inv_y)
inv_y = inv_y[:,0]

y_actual = pd.DataFrame(inv_y, columns=['Bearing Pressure (kPa)'])
y_hat = pd.DataFrame(inv_yhat, columns=['Predicted Bearing Pressure (kPa)'])
plt.figure(figsize=(11, 6))
plt.plot(y_actual, linestyle='solid', color='r')
plt.plot(y_hat, linestyle='dashed', color='b')
plt.legend(['Actual','Predicted'], loc='best', prop={'size': 14})
plt.title('Bearing Pressure ', weight='bold', fontsize=16)
plt.ylabel('(kPa)', weight='bold', fontsize=14)
plt.xlabel('Minutes Set no.', weight='bold', fontsize=14)
plt.xticks(weight='bold', fontsize=12, rotation=45)
plt.yticks(weight='bold', fontsize=12)
plt.grid(color = 'y', linewidth='0.5')
plt.show()

# calculate RMSE
rmse = sqrt(mean_squared_error(inv_y, inv_yhat))
print('Test RMSE: %.3f' % rmse)

Test RMSE: 30.454

errors = inv_y-inv_yhat

rmse = (errors**2).mean()**0.5
rmse

30.454036912334175

import seaborn as sns

sns.distplot(errors, bins=5, kde=True);

Section 2 PowerDraw

Due to the discussion of the correlations, we requiere simply develop a model to predict the Power Draw from the Speed input parameter

The librararies for the section 2

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

import seaborn as sns
import scipy.stats as stats

df_train = pd.read_csv("data/sag_data_train.csv", index_col="Time", parse_dates=True)
df = df_train.iloc[:, 1:3]
df = df.reset_index(drop=True)

df.head()

	Power Draw (MW)	Speed (RPM)
0	11.737357	7.843532
1	11.543755	7.681607
2	11.169525	7.514173
3	11.035091	7.592248
4	11.418827	7.784895

dataset=df

dataset.plot(x='Speed (RPM)', y='Power Draw (MW)', style='o')
plt.title('Speed (RPM) vs Power Draw (MW)')
plt.xlabel('Speed (RPM)')
plt.ylabel('Power Draw (MW)')
plt.show()

It is interesting it seems that there are several outliers that makes us a good model

So I define a function that removes those outliers

I will use interquartile range (IQR), to filter the outiers. In descriptive statistics, the interquartile range (IQR), also called the midspread, middle 50%, or H‑spread, is a measure of statistical dispersion, being equal to the difference between 75th and 25th percentiles, or between upper and lower quartiles.IQR = Q3 − Q1. In other words, the IQR is the first quartile subtracted from the third quartile; these quartiles can be clearly seen on a box plot on the data. It is a trimmed estimator, defined as the 25% trimmed range, and is a commonly used robust measure of scale. Source: Wikipedia

def remove_outliers_irq(data,parameterX):
    import seaborn as sns
    Q1 = data[parameterX].quantile(0.25)
    Q3 = data[parameterX].quantile(0.75)
    IQR = Q3 - Q1    #IQR is interquartile range.
    filter = (data[parameterX] >= Q1 - 1.5 * IQR) & (data[parameterX] <= Q3 + 1.5 *IQR)
    newdf=data.loc[filter]  
    return newdf

Let see them the outlier in another plot

sns.boxplot(x=dataset['Power Draw (MW)'])

<matplotlib.axes._subplots.AxesSubplot at 0x2afdfd34ac0>

sns.boxplot(x=dataset['Speed (RPM)'])

<matplotlib.axes._subplots.AxesSubplot at 0x2afe179ef10>

data2=remove_outliers_irq(dataset,'Power Draw (MW)')

data3=remove_outliers_irq(data2,'Speed (RPM)')

sns.boxplot(x=data3['Power Draw (MW)'])

<matplotlib.axes._subplots.AxesSubplot at 0x2afe17f14f0>

sns.boxplot(x=data3['Speed (RPM)'])

<matplotlib.axes._subplots.AxesSubplot at 0x2afe1cbef10>

data3.plot(x='Speed (RPM)', y='Power Draw (MW)', style='o')
plt.title('Speed (RPM) vs Power Draw (MW)')
plt.xlabel('Speed (RPM)')
plt.ylabel('Power Draw (MW)')
plt.show()

So now the data its a bit better to use develop a model

I would like use the first column as X and Y the second what I want to predict

neworder = ['Speed (RPM)','Power Draw (MW)']
data3=data3.reindex(columns=neworder)

dataset=data3

We check if our data its okay to use, ( to use regression we have to know that we dont have any NAN or strings, we need numeric data)

dataset.isnull().any()

Speed (RPM)        False
Power Draw (MW)    False
dtype: bool

So its OKAY

dataset.head()

	Speed (RPM)	Power Draw (MW)
0	7.843532	11.737357
1	7.681607	11.543755
2	7.514173	11.169525
3	7.592248	11.035091
4	7.784895	11.418827

FITTING PROCEDURE

I need take the numbers from the dataset

X = dataset.iloc[:, :-1].values
y = dataset.iloc[:, 1].values

We use the sklearn to create the model of LinearRegression

from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)

from sklearn.linear_model import LinearRegression
regressor = LinearRegression()
regressor.fit(X_train, y_train)

LinearRegression()

Finally I got the the model

y_pred = regressor.predict(X_test)

from sklearn import metrics
print('Mean Absolute Error:', metrics.mean_absolute_error(y_test, y_pred))
print('Mean Squared Error:', metrics.mean_squared_error(y_test, y_pred))
print('Root Mean Squared Error:', np.sqrt(metrics.mean_squared_error(y_test, y_pred)))

Mean Absolute Error: 0.36206541600380093
Mean Squared Error: 0.2078778493453066
Root Mean Squared Error: 0.45593623385875637

Prediction

I repeat the same procedure used for the first section

df_test = pd.read_csv("data/sag_data_test.csv", index_col="Time", parse_dates=True)
df_test_5min=df_test['2015-12-21 00:00:00':'2015-12-21 00:05:00']
dframe=df_test_5min.iloc[1:]
actual_power=dframe['Power Draw (MW)'].values # this  has shape (XXX, ) - It's 1D
input_Speed=dframe['Speed (RPM)'].values.reshape(-1,1) # this  has shape (XXX, 1) - it's 2D

I perform the prediction

PowerDraw_prediction = regressor.predict(input_Speed)
y_hat2 = pd.DataFrame(PowerDraw_prediction, columns=['Predicted Power Draw (MW)']) # I create the dataframe of the prediction

plt.figure(figsize=(11, 6))
plt.plot(actual_power, linestyle='solid', color='r')
plt.plot(y_hat2, linestyle='dashed', color='b')
plt.legend(['Actual','Predicted'], loc='best', prop={'size': 14})
plt.title('Predicted Power Draw  ', weight='bold', fontsize=16)
plt.ylabel('(MW)', weight='bold', fontsize=14)
plt.xlabel('Minutes Set no.', weight='bold', fontsize=14)
plt.xticks(weight='bold', fontsize=12, rotation=45)
plt.yticks(weight='bold', fontsize=12)
plt.grid(color = 'y', linewidth='0.5')
plt.show()

from math import sqrt
from sklearn.metrics import mean_squared_error
# calculate RMSE
rmse = sqrt(mean_squared_error(actual_power, y_hat2))
print('Test RMSE: %.3f' % rmse)

Test RMSE: 0.239

pred_power=y_hat2['Predicted Power Draw (MW)'].values
errors2 = actual_power -  pred_power
rmse = (errors2**2).mean()**0.5
rmse

0.23905494524804163

sns.distplot(errors2, bins=6, kde=True);

Conclusion

We have developed two models, one for the prediction of Bearing Pressure and another Power Draw For the Bearing Pressure we have used as a model a neural network LSTM and for Power Draw LinearRegresor

For the first five minutes we create a dataframe called result

dframe=df_test_5min.iloc[1:]

result = pd.concat([ y_hat,y_hat2,dframe.reset_index()], axis=1, sort=False)

result.set_index('Time', inplace=True)

result

	Predicted Bearing Pressure (kPa)	Predicted Power Draw (MW)	Bearing Pressure (kPa)	Power Draw (MW)	Speed (RPM)	Dilution Flow Rate (m3/h)	Conveyor Belt Feed Rate (t/h)	Conveyor Belt PSD Fines (%)
Time
2015-12-21 00:01:00	5363.178223	12.126732	5343.775163	12.491533	7.889129	1057.823776	3042.468414	44.751640
2015-12-21 00:02:00	5354.933594	12.314203	5350.858199	12.290888	7.970280	1036.836364	2959.495635	42.130604
2015-12-21 00:03:00	5351.686035	13.049788	5361.138694	12.786268	8.288696	1054.944698	2994.353878	38.415417
2015-12-21 00:04:00	5355.383789	13.033396	5294.834751	13.010269	8.281600	1084.639678	2951.411343	45.020363
2015-12-21 00:05:00	5333.299805	12.544739	5311.194279	12.831332	8.070073	1078.903247	2995.219207	40.796616

result .plot(subplots=True, layout=(-1, 1), lw=1, figsize=(12, 6))
plt.tight_layout()

Congratulations! we could develop a forecast by using neural networks LSTM .