composite MLP models using multiprocessing

The model can be solved as an exponential function. It can be modelled by a neural network.

Using the following network parameters, a [coefficient of determination][2] of 0.953 was obtained on the data.

The following network parameters are used:

* Inputs: `TempCondOut`, `PerFreqFanA`, `PerFreqFanB`, `TempAmbient`, `TempWetBulb`

* Output: `TempCondIn`

momentum: 0.9

```

Fan power control signals are bi-modally distributed with low variance arount 100% and 0% (see [trends][3]). A minority of control signals fall in the (0%, 95%) interval. This may cause the model to simply learn system dynamics for the modes of the distribution. Two approaches are used:

#### Single MLP Model

A single MLP is trained on the entirety of the data. This achieves a [coefficient of determination][2] of 0.953.

#### Composite MLP

Three identical MLPs are trained separately on clusters of samples where the control signals, `PerFreqFanA` and `PerFreqFanB` are:

* Equal to 0

* Between 0 and 0.95

* Greater than 0.95

The following results are obtained:

| Cluster      	| Coefficient of determination 	|

|--------------	|------------------------------	|

| == 0         	| 0.97                         	|

| < 0 & < 0.95 	| 0.67                         	|

| > 0.95       	| 0.87                         	|

Giving a weighed coefficient of determination of 0.86.

[1]: 0-thermo-basics.md

[2]: https://en.wikipedia.org/wiki/Coefficient_of_determination

[3]: 6-trends.md

 No newline at end of file

%% Cell type:code id: tags:

``` python

%matplotlib notebook

%reload_ext autoreload

%autoreload 2

import datetime

from os import path, environ

from copy import deepcopy

from multiprocessing import Pool

import numpy as np

import pandas as pd

from sklearn.preprocessing import StandardScaler

from sklearn.pipeline import Pipeline

from sklearn.neural_network import MLPRegressor

from sklearn.model_selection import train_test_split

import ipyvolume as ipv

from models import fit_composite_model

# source file, see docs/5-dataset.md for info on field names

chiller_file = path.join(environ['DATADIR'], 'EngineeringScienceBuilding', 'Chillers.csv')

plot_path = path.join('..', 'docs', 'img')

```

%% Cell type:code id: tags:

``` python

# Read pre-processed data

df = pd.read_csv(chiller_file, index_col='Time', parse_dates=True, dtype=float)

df.dropna(inplace=True)

```

%% Cell type:markdown id: tags:

# Evaporative Cooling Model

Modelling the relationship between the cooling done by the cooling tower and the environmental, system, and control inputs.

$$

\begin{align*}

T(t) &= T(0) e^{-\frac{k T_a v_f R}{T_w c_m m} t}

\end{align*}

$$

Where:

* $T(0)$ is the warm entering water temperature `TempCondOut`,

* $T(t)$ is the exiting cool water temperature `TempCondIn`,

* $T_a$ is ambient air temperature `TempAmbient`,

* $v_f$ is avarage fan speed `0.5 * (PerFreqFanA + PerFreqFanB)`,

* $R$ is solar irradiance - constant if assuming shade.

* $T_w$ is wet-bulb temperature `TempWetBulb`,

* $c_m$ is specific mass heat capacity of water,

* $m$ is the mass of water being cooled - constant if assuming steady flow rate.

The model predicts $T(t)$ from all other factors.

## Multi-Layer Perceptron

Assuming steady flow rate, constant irradiance.

The motivation for an MLP model is the assumption that inputs manifest instantaneously as outputs with not delay. This is a simplification as water takes some time to cycle through the cooling tower.

The inputs are chosen assuming steady flow rate, constant irradiance.

The control variables are `PerFreqFanA` and `PerFreqFanB` which usually track each other. They are distributed bi-modally around 100% and 0% power with a minority of samples falling in between. A concern is that the model may learn to treat control inputs as constant. Two approaches are chosen:

* A single model is learned over all data,

* Samples are clustered by control=0, control >= 0.95, and 0 < control  < 0.95. A separate model is learned for each cluster.

%% Cell type:code id: tags:

``` python

# Data

feature_cols = ['TempCondOut', 'PerFreqFanA', 'PerFreqFanB', 'TempAmbient', 'TempWetBulb']

X = df.loc[:, feature_cols]

Y = df['TempCondIn']

X, Y = df.loc[:, feature_cols], df['TempCondIn']

# setting up cluster selectors

c1 = df['PerFreqFanA'] >= 0.95

c0 = df['PerFreqFanA'] == 0

cmid = ~ (c1 | c0)

# generating clusters

X1, Y1 = X[c1], Y[c1]

Xmid, Ymid = X[cmid], Y[cmid]

X0, Y0 = X[c0], Y[c0]

# generating training/testing sets for each cluster

Xtrain, Xtest, Ytrain, Ytest = train_test_split(X, Y, test_size=0.1)

X1train, X1test, Y1train, Y1test = train_test_split(X1, Y1, test_size=0.1)

Xmidtrain, Xmidtest, Ymidtrain, Ymidtest = train_test_split(Xmid, Ymid, test_size=0.1)

X0train, X0test, Y0train, Y0test = train_test_split(X0, Y0, test_size=0.1)

```

%% Cell type:markdown id: tags:

### Single MLP

%% Cell type:code id: tags:

``` python

# Pipeline

scaler = StandardScaler()              # scale to 0 mean and unit variance

regressor = MLPRegressor(hidden_layer_sizes=(20,20),

                         learning_rate_init=1e-3,

                         verbose=True) # regression model

reg = Pipeline([('scaler', scaler), ('regressor', regressor)])

est = Pipeline([('scaler', scaler), ('regressor', regressor)])

```

%% Cell type:code id: tags:

``` python

est.fit(Xtrain, Ytrain)

```

%% Cell type:code id: tags:

``` python

print('Score on all control:\t{:.4f}'.format(est.score(Xtest, Ytest)))

print('Score on >95%-control:\t{:.4f}'.format(est.score(X0test, Y0test)))

print('Score on mid-control:\t{:.4f}'.format(est.score(X1test, Y1test)))

print('Score on 0%-control:\t{:.4f}'.format(est.score(Xmidtest, Ymidtest)))

```

%% Cell type:markdown id: tags:

### Composite MLP

%% Cell type:code id: tags:

``` python

reg.fit(Xtrain, Ytrain)

reg.score(Xtest, Ytest)

scaler = StandardScaler()              # scale to 0 mean and unit variance

regressor = MLPRegressor(hidden_layer_sizes=(20,20),

                         learning_rate_init=1e-3,

                         verbose=False,

                         max_iter=1000,

                         solver='adam') # regression model

est = Pipeline([('scaler', scaler), ('regressor', regressor)])

est1, estmid, est0 = fit_composite_model(est, zip((X1train, Xmidtrain, X0train),

                                                  (Y1train, Ymidtrain, Y0train)))

print('Loss on >95%-control:\t{:.4f}'.format(est1.named_steps['regressor'].loss_))

print('Loss on mid-control:\t{:.4f}'.format(estmid.named_steps['regressor'].loss_))

print('Loss on 0%-control:\t{:.4f}'.format(est0.named_steps['regressor'].loss_))

```

%% Cell type:code id: tags:

``` python

# Plotting learned function

score1 = est1.score(X1test, Y1test)

scoremid = estmid.score(Xmidtest, Ymidtest)

score0 = est0.score(X0test, Y0test)

scoreall = (len(X1test)*score1 + len(Xmidtest)*scoremid + len(X0test)*score0) / len(Xtest)

print('Score on all control:\t{:.4f}'.format(scoreall))

print('Score on >95%-control:\t{:.4f}'.format(score1))

print('Score on mid-control:\t{:.4f}'.format(scoremid))

print('Score on 0%-control:\t{:.4f}'.format(score0))

```

%% Output

    Score on all control:	0.8641

    Score on >95%-control:	0.8700

    Score on mid-control:	0.6506

    Score on 0%-control:	0.9793

%% Cell type:code id: tags:

``` python

%matplotlib inline

%reload_ext autoreload

%autoreload 2

import datetime

from os import path, environ

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

from sklearn.preprocessing import MinMaxScaler, StandardScaler

from sklearn.pipeline import Pipeline

from sklearn.linear_model import LinearRegression

from sklearn.model_selection import train_test_split

from sklearn.cluster import DBSCAN

from sklearn.metrics import mutual_info_score

import ipyvolume as ipv

from thermo import CONSTANTS

from preprocess import POW_FIELDS

# source file, see docs/5-dataset.md for info on field names

Chiller1File = path.join(environ['DATADIR'], 'EngineeringScienceBuilding', 'Chillers.csv')

plot_path = path.join('..', 'docs', 'img')

```

%% Cell type:code id: tags:

``` python

# Read pre-processed data

df = pd.read_csv(Chiller1File, index_col='Time', parse_dates=True, dtype=float)

df.dropna(inplace=True)

```

%% Cell type:markdown id: tags:

# Correlations

## Cross-correlation

Pearson product-moment correlation coefficients between variables. Measure of linear relationship between variables.

$$

r_{X,Y} = \frac{\textrm{cov}(X,Y)}{\sigma_X \sigma_Y}

$$

%% Cell type:code id: tags:

``` python

corr = np.corrcoef(df, rowvar=False)

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

axImg = ax.matshow(corr, vmin=-1, vmax=1)

fig.colorbar(axImg)

plt.xticks(np.arange(len(df.columns)), df.columns, rotation='vertical')

plt.yticks(np.arange(len(df.columns)), df.columns);

plt.savefig(path.join(plot_path, '7-pearson-correlation.png'))

```

%% Cell type:markdown id: tags:

### Mutual information

## Mutual information

Mutual information is a measure of how much information about the distribution variable $X$ is contained in the distribution of variable $Y$.

$$

\mathrm{MI}(X,Y) = \sum_{i}^{\mid X \mid} \sum_{j}^{\mid Y \mid} \frac{\mid X_i \cap Y_j \mid}{N} \log{\frac{N \; \mid X_i \cap Y_j \mid}{\mid X_i \mid \; \mid Y_j \mid}}

$$

Where $N$ is the total number of samples, $i$ and $j$ are class labels for samples (in this case, histogram bins). $X_i$ is the number of samples of $X$ with a label $i$. The MI score is normalized to lie between 0 and 1.

%% Cell type:code id: tags:

``` python

bins = 32

mi = np.empty((len(df.columns), len(df.columns)), dtype=float)

# Standardizing values for 0 mean and unit variance

scaled = StandardScaler().fit_transform(df)

# Calculating MI matrix for upper triangular half explicitly,

# the lower half is symmetric.

for i, col1 in enumerate(df.columns):     # row

    for j, col2 in enumerate(df.columns): # column

        if j < i:

            continue

        else:

            hist,_ , _ = np.histogram2d(scaled[:,i], scaled[:,j], bins=bins, density=False)

            N = hist.sum()

            # Normalize MI, taken from source of sklearn's

            # normalized_mutual_info_score to work with a

            # contingency table.

            p_i, p_j = np.sum(hist, axis=0) / N, np.sum(hist, axis=1) / N

            non_zero_i, non_zero_j = p_i > 0, p_j > 0

            ent_i = -np.sum(p_i[non_zero_i] * np.log(p_i[non_zero_i]))

            ent_j = -np.sum(p_j[non_zero_j] * np.log(p_j[non_zero_j]))

            norm = 0.5 * (ent_i + ent_j)

            mi[i,j] = mutual_info_score(None, None, contingency=hist) / norm

            mi[j,i] = mi[i,j]

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

axImg = ax.matshow(mi)

fig.colorbar(axImg)

plt.xticks(np.arange(len(df.columns)), df.columns, rotation='vertical')

plt.yticks(np.arange(len(df.columns)), df.columns);

plt.savefig(path.join(plot_path, '7-mutual-information.png'))

```

%% Cell type:markdown id: tags:

## Fan speed & fan power

## Temporal correlation

*Hypothesis*: Fan power depends on ambient temperature `TempAmbient`, relative humidity `PerHumidity`, and fan speed setting `PerFreqFan[A | B]`.

*Result*: Very little relationship from current data. Primarily because fan speed is ~100%. Need more variation.

%% Cell type:code id: tags:

``` python

# Data

X = df.loc[:, ('TempAmbient', 'PerHumidity', 'PerFreqFanA')]

Y = df.loc[:, 'PowFanA']

Xtrain, Xtest, Ytrain, Ytest = train_test_split(X, Y, test_size=0.1)

```

%% Cell type:code id: tags:

``` python

# Pipeline

scaler = MinMaxScaler()        # scale to [0-1] range

regressor = LinearRegression() # regression model

reg = Pipeline([('scaler', scaler), ('regressor', regressor)])

```

%% Cell type:code id: tags:

``` python

reg.fit(Xtrain, Ytrain)

reg.score(Xtest, Ytest)

```

%% Output

    0.1107982030593837

Measure the relationship variables after introducing a lag.

%% Cell type:code id: tags:

``` python

f = plt.figure(figsize=(8,8))

ax = f.add_subplot(111, projection='3d')

ax.scatter3D(X.iloc[:, 0], X.iloc[:, 1], Y, c=Y)

ax.set_xlabel('Ambient Temperature (K)')

ax.set_ylabel('Relative Humidity')

ax.set_zlabel('Fan Power (W)')

plt.savefig(path.join(plot_path, '7-fan-power-vs-temp-humidity.png'))

# For each lag:

#   For each variable:

#     lag all other variables

#

```

%% Cell type:markdown id: tags:

# Clusters

## Temperature

%% Cell type:code id: tags:

``` python

X = df.loc[:, ('TempAmbient', 'TempWetBulb')]

X['DeltaTemp'] = df['TempCondOut'] - df['TempCondIn']

labels = DBSCAN().fit_predict(X)

```

%% Cell type:code id: tags:

``` python

ipv.clear()

cmap = plt.cm.Accent(labels)

f = ipv.scatter(X['TempAmbient'], X['DeltaTemp'], X['TempWetBulb'], color=cmap,

            size=10, marker='point_2d')

ipv.xyzlim(280, 310)

ipv.ylim(0, 10)

ipv.xyzlabel('TempAmbient', 'DeltaTemp', 'TempWetBulb')

ipv.show()

```

%% Cell type:code id: tags:

``` python

def view(fig, n, frac):

    ipv.view(360*frac)

ipv.movie(path.join(plot_path, '7-ct-temp-clusters.gif'), view, fps=12, frames=36)

```

%% Output

%% Cell type:markdown id: tags:

![](../docs/img/7-ct-temp-clusers.gif)

%% Cell type:markdown id: tags:

## Power

%% Cell type:code id: tags:

``` python

X = df.loc[:, ('PowConP', 'PowFanA', 'PowFanB')]

labels = DBSCAN().fit_predict(X)

```

%% Cell type:code id: tags:

``` python

ipv.clear()

cmap = plt.cm.Accent(labels)

f = ipv.scatter(X['PowConP'], X['PowFanA'], X['PowFanB'], color=cmap,

            size=10, marker='point_2d')

# ipv.xyzlim(280, 310)

ipv.xyzlabel('PowConP', 'PowFanA', 'PowFanB')

ipv.show()

```

%% Cell type:code id: tags:

``` python

def view(fig, n, frac):

    ipv.view(360*frac)

ipv.movie(path.join(plot_path, '../docs/img/7-ct-power-clusters.gif'), view, fps=12, frames=36)

```

"""

Utility definitions for Models.ipynb

"""

from multiprocessing import Pool

from os import cpu_count

from typing import Iterable, Tuple

import numpy as np

from sklearn import clone

from sklearn.neural_network import MLPRegressor

def _fit(*args):

    """Multi-processing payload function used by fit_composite_model"""

    est, (x, y) = args

    return est.fit(x, y)

def fit_composite_model(estimator: MLPRegressor,

    data: Iterable[Tuple[np.ndarray, np.ndarray]]) -> Iterable[MLPRegressor]:

    """

    Fits copies of an estimator to different datasets in parallel.

    Args:

    * `estimator`: An estimator instance with a `fit(X, y)` method which returns

    the instance.

    * `data`: An iterable of tuples of arrays: [(train1, train2,..), (test1, test2,..)].

    Returns:

    * A list of fitted estimators.

    """

    data = list(data)

    estimators = [clone(estimator) for _ in data]

    with Pool(min(len(data), cpu_count())) as pool:

        return pool.starmap(_fit, zip(estimators, data))