Commit 719ac105 authored by hazrmard's avatar hazrmard
Browse files

bugfix with v1 preprocessing ambient temperatures, started RL control notebook

parent 7922b4c2

src/ICAART.ipynb→src/ICAART-2020.ipynb

+10 −20
Original line number Diff line number Diff line 
%% Cell type:code id: tags:



``` python

import sys, os

pytorchbridge_path = os.path.abspath('../../pyTorchBridge')

if pytorchbridge_path not in sys.path:

    sys.path.append(pytorchbridge_path)

```



%% Cell type:code id: tags:



``` python

%matplotlib notebook

%reload_ext autoreload

%autoreload 2



import datetime

from os import path, environ

import pickle



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, GridSearchCV, KFold

import torch

import torch.nn as nn

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

from pytorchbridge import TorchEstimator



from utils.stats import temporal_correlations

from plotting import model_surface, plot_surface

from controller import GridSearchController, \

                       BinaryApproachController, \

                       QuasiNewtonController

# 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

# Data selection 'all' or 'chiller_on' or 'fan_on'

MODE = 'chiller_on'

# Read pre-processed data:

# Pytorch uses float32 as default type for weights etc,

# so input data points are also read in the same type.

df = pd.read_csv(chiller_file, index_col='Time',

                 parse_dates=['Time'], dtype=np.float32)

df.dropna(inplace=True)

if MODE == 'chiller_on':

    df = df[df['PowChi'] != 0.]

if MODE == 'fan_on':

    df = df[(df['PerFreqFanA'] != 0.) | df['PerFreqFanB'] != 0.]

print(len(df), 'Records')

```



%% Cell type:code id: tags:



``` python

df.TempCondOut.plot.hist()

```



%% Cell type:markdown id: tags:



# 1. Post-chiller temp model



%% Cell type:code id: tags:



``` python

# Data

feature_cols = ('TempCondIn',                 'TempEvapIn',

                'TempEvapOut', 'TempAmbient', 'TempWetBulb',

                'FlowEvap', 'PowConP')

target_cols = ('TempCondOut',)

# Normalizing data to have 0 mean and 1 variance

XChi, YChi = df.loc[:, feature_cols], df.loc[:,target_cols]

ScalerXChi, ScalerYChi = StandardScaler().fit(XChi), StandardScaler().fit(YChi)

XChi, YChi = ScalerXChi.transform(XChi), np.squeeze(ScalerYChi.transform(YChi))

XChitrain, XChitest, YChitrain, YChitest = train_test_split(XChi, YChi, test_size=0.1)



norm_chi_var = lambda i, v: ScalerXChi.mean_[i] + np.sqrt(ScalerXChi.var_[i]) * v

norm_chi = lambda v: ScalerYChi.inverse_transform(v)

std_chi_var = lambda i, v: (v - ScalerXChi.mean_[i]) / np.sqrt(ScalerXChi.var_[i])



CHILLER_MODELS = {}

```



%% Cell type:markdown id: tags:



## 1.1 Multi-layer perceptron



%% Cell type:code id: tags:



``` python

# Searching parameter grid for best hyperparameters

param_grid = {

#     'hidden_layer_sizes': [(4,4), (8,8), (16,16), (8,8,8), (16,16,16)],

    'hidden_layer_sizes': [(16,16), (16,16,16)],

#     'learning_rate_init': [1e-2, 1e-3, 1e-4],

    'learning_rate_init': [1e-3, 1e-4],

    'hidden_layer_sizes': [(4,4), (8,8), (16,16), (8,8,8), (16,16,16)],

    'learning_rate_init': [1e-2, 1e-3, 1e-4],

    'activation': ['relu', 'tanh', 'logistic']

}



grid_search = GridSearchCV(MLPRegressor(), param_grid, n_jobs=4, verbose=1,

                           cv=KFold(3, shuffle=True))

grid_search.fit(XChitrain, YChitrain)



est = grid_search.best_estimator_

print('Test score:', est.score(XChitest, YChitest))



CHILLER_MODELS['SingleMLP'] = est

pickle.dump(est, open('./bin/chiller_{}.pickle'.format(MODE), 'wb'))



grid_res = pd.DataFrame(grid_search.cv_results_).sort_values('rank_test_score')

grid_res.head()

```



%% Cell type:code id: tags:



``` python

# Load model from file instead of training

est = pickle.load(open('./bin/chiller_{}.pickle'.format(MODE), 'rb'))

CHILLER_MODELS['SingleMLP'] = est

```



%% Cell type:code id: tags:



``` python

# Plot model predictions when varying 2 fields at a single time instant

dt = datetime.datetime(2018, 6, 2, 15, 0, 0)

select = df.index == dt

singleX, singleY = XChi[select], YChi[select]

var, var_idx = ('TempCondIn', 'PowConP'), (0, 6)

vary_range = ((std_chi_var(var_idx[0], 293), 1.),

              (std_chi_var(var_idx[1], 0), 1.5))



x, y, z = model_surface(CHILLER_MODELS['SingleMLP'], singleX, var_idx,

                        vary_range, (10, 10))

ax = plot_surface(norm_chi_var(var_idx[0], x),

                  norm_chi_var(var_idx[1], y),

                  norm_chi(z),

                  alpha=0.75, cmap=plt.get_cmap('coolwarm'))

ax.set_xlabel(var[0], labelpad=20)

ax.set_ylabel(var[1])

ax.set_zlabel('Temperature / K')

plt.title('Temperature Extrapolation ' + dt.isoformat());

```



%% Cell type:markdown id: tags:



# 2. Energy Model



%% Cell type:code id: tags:



``` python

# Data

feature_cols = ('TempCondIn', 'TempCondOut', 'TempEvapIn',

                'TempEvapOut', 'TempAmbient', 'TempWetBulb',

                'FlowEvap', 'PowConP')

target_cols = ('PowChi',)

# Normalizing data to have 0 mean and 1 variance

XEnergy, YEnergy = df.loc[:, feature_cols], df.loc[:,target_cols]

ScalerXEnergy, ScalerYEnergy = StandardScaler().fit(XEnergy), StandardScaler().fit(YEnergy)

XEnergy, YEnergy = ScalerXEnergy.transform(XEnergy), np.squeeze(ScalerYEnergy.transform(YEnergy))

XEnergytrain, XEnergytest, YEnergytrain, YEnergytest = train_test_split(XEnergy, YEnergy, test_size=0.1)



norm_energy_var = lambda i, v: ScalerXEnergy.mean_[i] + np.sqrt(ScalerXEnergy.var_[i]) * v

norm_energy = lambda v: ScalerYEnergy.inverse_transform(v)

std_energy_var = lambda i, v: (v - ScalerXEnergy.mean_[i]) / np.sqrt(ScalerXEnergy.var_[i])



ENERGY_MODELS = {}

```



%% Cell type:code id: tags:



``` python

# Measuring temporal correlations between features and targets`

lags = tuple(range(10))

corrs = temporal_correlations(XEnergy, lags=lags, targets=YEnergy)

maxcorridx = np.argmax(np.square(corrs), axis=0)

maxcorrs = np.asarray(lags)[maxcorridx.flatten()]

print('Feature vs most correlated lag')

for feature, lag in zip(feature_cols, maxcorrs):

    print('{:12s}:\t{}'.format(feature, lag))

```



%% Cell type:markdown id: tags:



## 2.1 Multi-layer perceptron



%% Cell type:code id: tags:



``` python

# Searching parameter grid for best hyperparameters

param_grid = {

#     'hidden_layer_sizes': [(4,4), (8,8), (16,16), (8,8,8), (16,16,16)],

    'hidden_layer_sizes': [(16,16), (16,16,16)],

#     'learning_rate_init': [1e-2, 1e-3, 1e-4],

    'learning_rate_init': [1e-3, 1e-4],

    'activation': ['relu', 'tanh', 'logistic']

}



grid_search = GridSearchCV(MLPRegressor(), param_grid, n_jobs=4,

                           verbose=1, cv=KFold(3, shuffle=True))

grid_search.fit(XEnergytrain, YEnergytrain)



est = grid_search.best_estimator_

print('Test score:', est.score(XEnergytest, YEnergytest))



ENERGY_MODELS['SingleMLP'] = est

pickle.dump(est, open('./bin/energy_{}.pickle'.format(MODE), 'wb'))



grid_res = pd.DataFrame(grid_search.cv_results_)\

             .sort_values('rank_test_score')

grid_res.head()

```



%% Cell type:code id: tags:



``` python

# Load model from file instead of training

est = pickle.load(open('./bin/energy_{}.pickle'.format(MODE), 'rb'))

ENERGY_MODELS['SingleMLP'] = est

```



%% Cell type:code id: tags:



``` python

# Plot prediction vs actual across day or timestamp

# time = datetime.time(12,0,0)

time = datetime.date(2018,6,2)



if isinstance(time, datetime.time):

    select = df.index.time == time

    tseries = df.index.date[select]

elif isinstance(time, datetime.date):

    select = df.index.date == time

    tseries = df.index.time[select]

Xfiltered = XEnergy[select]

Yfiltered = YEnergy[select]



plt.scatter(np.arange(len(tseries)),

            ScalerYEnergy.inverse_transform(Yfiltered),

            c='r', label='Historic')

plt.scatter(np.arange(len(tseries)),

            ScalerYEnergy.inverse_transform(est.predict(Xfiltered)),

            c='g', label='Predicted')

plt.title('Power Model ' + time.isoformat())

plt.ylabel('Power / W')

plt.xlabel('Date/Time')

plt.legend();

ticks, labels = plt.xticks()

ticks = np.asarray([t for t in ticks if 0 <= t < len(tseries)])

plt.xticks(ticks, tseries[ticks.astype(int)]);

```



%% Cell type:code id: tags:



``` python

# Plot model predictions when varying a field

# var= 'TempCondIn'

var= 'PowConP'

var_idx = feature_cols.index(var)



vary_range = (std_energy_var(var_idx, 273), 2.)

_, y, z = model_surface(ENERGY_MODELS['SingleMLP'], Xfiltered, (var_idx,),

                        (vary_range,), (20,))

ax = plot_surface(tseries,

                  norm_energy_var(var_idx, y),

                  norm_energy(z),

                  alpha=0.75, cmap=plt.get_cmap('coolwarm'))

ax.set_xlabel('Date/Time', labelpad=20)

ax.set_ylabel(var)

ax.set_zlabel('Power / W')

plt.title('Power Extrapolation ' + time.isoformat());

```



%% Cell type:code id: tags:



``` python

# Plot model predictions when varying 2 fields at a single time instant

dt = datetime.datetime(2018, 6, 2, 15, 0, 0)

select = df.index == dt

singleX, singleY = XEnergy[select], YEnergy[select]

var, var_idx = ('TempCondIn', 'PowConP'), (0, 7)

vary_range = ((std_energy_var(var_idx[0], 293), 1.),

              (std_energy_var(var_idx[1], 0), 1.))



x, y, z = model_surface(ENERGY_MODELS['SingleMLP'], singleX, var_idx,

                        vary_range, (10, 10))

ax = plot_surface(norm_energy_var(var_idx[0], x),

                  norm_energy_var(var_idx[1], y),

                  norm_energy(z),

                  alpha=0.75, cmap=plt.get_cmap('coolwarm'))

ax.set_xlabel(var[0], labelpad=20)

ax.set_ylabel(var[1])

ax.set_zlabel('Power / W')

plt.title('Power Extrapolation ' + dt.isoformat());

```



%% Cell type:markdown id: tags:



# 3. Evaporative Cooling Model



%% Cell type:code id: tags:



``` python

# Data

fan_cols = ('PerFreqFanA', 'PerFreqFanB')

feature_cols = ('TempCondOut', 'TempAmbient', 'TempWetBulb', 'PowConP')

target_cols = ('TempCondIn',)

# Combining 2 fan speed controls into a single variable (averaged)

X = pd.concat((df.loc[:, fan_cols].mean(axis=1),

               df.loc[:, feature_cols]), axis=1)

Y = df.loc[:,target_cols]

# Normalizing data to have 0 mean and 1 variance

ScalerX, ScalerY = StandardScaler().fit(X), StandardScaler().fit(Y)

X, Y = ScalerX.transform(X), np.squeeze(ScalerY.transform(Y))

# generating training/testing sets

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



# Convert feature variables or target to normal or standardized form:

# i - is the index of the variable in the feature array

# v - is an array of values to transform

norm_var = lambda i, v: ScalerX.mean_[i] + np.sqrt(ScalerX.var_[i]) * v

std_var = lambda i, v: (v - ScalerX.mean_[i]) / np.sqrt(ScalerX.var_[i])

norm_temp = lambda v: ScalerY.inverse_transform([v])[0]

std_temp = lambda v: ScalerY.transform([v])[0]



EVAP_MODELS = {}

```



%% Cell type:code id: tags:



``` python

# Distribution of control variables

plt.figure(figsize=(10,4))

plt.subplot(1,2,1)

df['PowConP'].hist(bins=20);

plt.xlabel('PowConP / W')

plt.ylabel('Frequency')

plt.ylim(0, 60000)

plt.yticks(rotation=70)

plt.subplot(1,2,2)

((df['PerFreqFanA'] + df['PerFreqFanB']) / 2).hist(bins=20)

plt.xlabel('PerFreqFan(A | B)')

plt.ylim(0, 60000)

plt.ylabel('Frequency')

plt.yticks(rotation=70);



# Fraction of scarce values

sparse_fan = (std_var(0, 0.1) < Xtrain[:, 0]) & (Xtrain[:, 0] < std_var(0, 0.9))

sparse_pump = (std_var(4, 0.1 * max(Xtrain[:, 4])) < Xtrain[:, 4]) & (Xtrain[:, 4] < std_var(0, 0.9 * max(Xtrain[:, 4])))

sparse_only_fan = sparse_fan & ~sparse_pump

sparse_only_pump = ~sparse_fan & sparse_pump

sparse_joint = sparse_fan & sparse_pump



print('Sparsity - fan:', sum(sparse_fan) / len(sparse_fan))

print('Sparsity - only fan:', sum(sparse_only_fan) / len(sparse_only_fan))

print('Sparsity - pump:', sum(sparse_pump) / len(sparse_pump))

print('Sparsity - only pump:', sum(sparse_only_pump) / len(sparse_only_pump))

print('Sparsity - joint:', sum(sparse_joint) / len(sparse_pump))

```



%% Cell type:code id: tags:



``` python

# Re-distributing variables for more uniform spread

sparsity = int(len(sparse_joint) / sum(sparse_joint))

Xtrain_sparse = np.tile(Xtrain[sparse_joint], (sparsity - 1, 1))

Ytrain_sparse = np.tile(Ytrain[sparse_joint], (sparsity - 1))

Xtrain = np.vstack((Xtrain, Xtrain_sparse))

Ytrain = np.hstack((Ytrain, Ytrain_sparse))

print(Xtrain.shape, Ytrain.shape)



# Distribution of control variables

plt.figure(figsize=(10,4))

plt.subplot(1,2,1)

pd.DataFrame({'PowConP': norm_var(4, Xtrain[:, 4])}).hist(bins=20, ax=plt.gca())

plt.xlabel('PowConP / W')

plt.ylabel('Frequency')

plt.title('')

plt.ylim(0, 60000)

plt.yticks(rotation=70)

plt.subplot(1,2,2)

pd.DataFrame({'PerFreqFan': norm_var(0, Xtrain[:, 0])}).hist(bins=20, ax=plt.gca())

plt.xlabel('PerFreqFan(A | B)')

plt.ylim(0, 60000)

plt.ylabel('Frequency')

plt.title('')

plt.yticks(rotation=70);

```



%% Cell type:code id: tags:



``` python

# Measuring temporal correlations between features and targets`

lags = tuple(range(20))

corrs = temporal_correlations(X, lags=lags, targets=Y)

maxcorridx = np.argmax(np.square(corrs), axis=0)

maxcorrs = np.asarray(lags)[maxcorridx.flatten()]

print('Feature vs most correlated lag')

for feature, lag in zip(('FanFreq',) + feature_cols, maxcorrs):

    print('{:12s}:\t{}'.format(feature, lag))

```



%% Cell type:markdown id: tags:



## 3.1 Multi-Layer Perceptron



%% Cell type:code id: tags:



``` python

# Searching parameter grid for best hyperparameters

param_grid = {

#     'hidden_layer_sizes': [(4, 4, 4), (8, 8, 8), (16, 16, 16)],

    'hidden_layer_sizes': [(16, 16), (16, 16, 16)],

#     'learning_rate_init': [1e-2, 1e-3, 1e-4],

    'learning_rate_init': [1e-3, 1e-4],

    'activation': ['relu', 'tanh', 'logistic']

}



grid_search = GridSearchCV(MLPRegressor(), param_grid, n_jobs=4, verbose=1, cv=KFold(3, shuffle=True))

grid_search.fit(Xtrain, Ytrain)



est = grid_search.best_estimator_

print('Test score:', est.score(Xtest, Ytest))



EVAP_MODELS['SingleMLP'] = est

pickle.dump(est, open('./bin/evap_{}.pickle'.format(MODE), 'wb'))



grid_res = pd.DataFrame(grid_search.cv_results_).sort_values('rank_test_score')

grid_res.head()

```



%% Cell type:code id: tags:



``` python

est = pickle.load(open('./bin/evap_{}.pickle'.format(MODE), 'rb'))

EVAP_MODELS['SingleMLP'] = est

```



%% Cell type:code id: tags:



``` python

# Get predictions

# time = datetime.time(4,0,0)

time = datetime.date(2018,6,2)



if isinstance(time, datetime.time):

    select = df.index.time == time

    tseries = df.index.date[select]

elif isinstance(time, datetime.date):

    select = df.index.date == time

    tseries = df.index.time[select]

Xfiltered = X[select]

Yfiltered = Y[select]

```



%% Cell type:code id: tags:



``` python

# Plot model predictions when varying a single field over a series of times

# var, var_idx = 'PerFreqFan', 0

var, var_idx = 'PowConP', 4



vary_range = (std_var(var_idx, 0), 1.)  # vary from 0 to 1 std dev above mean



_, y, z = model_surface(EVAP_MODELS['SingleMLP'], Xfiltered, (var_idx,), (vary_range,), (10,))

ax = plot_surface(tseries,

                  norm_var(var_idx, y),

                  norm_temp(z),

                  alpha=0.75, cmap=plt.get_cmap('coolwarm'))

ax.set_xlabel('Date/Time', labelpad=20)

ax.set_ylabel(var)

ax.set_zlabel('Temperature')

plt.title('Temperature Extrapolation ' + time.isoformat());

```



%% Cell type:code id: tags:



``` python

# Plot model predictions when varying 2 fields at a single time instant

dt = datetime.datetime(2018, 6, 2, 15, 0, 0)

select = df.index == dt

singleX, singleY = X[select], Y[select]

var, var_idx = ('PerFreqFan', 'PowConP'), (0, 4)

vary_range = ((std_var(var_idx[0], 0), 1.), (std_var(var_idx[1], 0), 1.))  # vary from 0 to 1 std dev above mean



x, y, z = model_surface(EVAP_MODELS['SingleMLP'], singleX, var_idx, vary_range, (10, 10))

ax = plot_surface(norm_var(var_idx[0], x),

                  norm_var(var_idx[1], y),

                  norm_temp(z),

                  alpha=0.75, cmap=plt.get_cmap('coolwarm'))

ax.set_xlabel(var[0], labelpad=20)

ax.set_ylabel(var[1])

ax.set_zlabel('TempCondIn')

plt.title('Temperature Extrapolation ' + dt.isoformat());

```



%% Cell type:code id: tags:



``` python

# Plot optimal control

vary_range = ((std_var(0, 0), 1.), (std_var(4, 2000), 1.))

# vary_range = ((std_var(0, 0), 1.),)

c = GridSearchController(est, vary_range, (10, 10), (0, 4))

# c = GridSearchController(est, vary_range, (10,), (0,))

# c = QuasiNewtonController(est, vary_range, (10, 10), (0, 4))

control = c.predict(Xfiltered)



plt.plot(control[:, 0], label='Optimal fan', ls=':', c='r')

plt.plot(Xfiltered[:, 0], label='Historic fan', ls=':', c='g')

plt.plot(control[:, 1], label='Optimal pump', ls='-', c='r')

plt.plot(Xfiltered[:, 4], label='Historic pump', ls='-', c='g')

plt.legend()

```



%% Cell type:markdown id: tags:



# 3. Optimal control



%% Cell type:code id: tags:



``` python

# time = datetime.time(14,0,0)

time = datetime.date(2018,6,2)

time = datetime.date(2018,4,2)

if isinstance(time, datetime.time):

    select = df.index.time == time

    tseries = df.index.date[select]

elif isinstance(time, datetime.date):

    select = df.index.date == time

    tseries = df.index.time[select]

Xfiltered = X[select]

Yfiltered = Y[select]

XEnergyfiltered = XEnergy[select]

YEnergyfiltered = YEnergy[select]



# Whether or not to use pump power as a control variable

USE_PUMP = True



if USE_PUMP:

    vary_range = ((std_var(0, 0), std_var(0, 1)),

                  (std_var(4, 2000), 1))

else:

    vary_range = ((std_var(0, 0), std_var(0, 1)),)





# The wetbulb + margin threshold temp used by the binary controller to

# switch fans on/off:

margin = (std_temp([3]) - std_temp([0]))[0]



est = EVAP_MODELS['SingleMLP']

est_energy = ENERGY_MODELS['SingleMLP']

est_chi = CHILLER_MODELS['SingleMLP']

if USE_PUMP:

    gridc = GridSearchController(est, vary_range, (10, 10), (0, 4))

else:

    gridc = GridSearchController(est, vary_range, (10,), (0,))

binaryc = BinaryApproachController(est, vary_range[0], 10, 0, margin)

```



%% Cell type:code id: tags:



``` python

# Comparison of control over a contiguous time period

Xdata = Xfiltered

Ydata = Yfiltered



baseline = np.squeeze(std_temp(norm_var(3, Xdata[:, 3])))



grid_control, binary_control = gridc.predict(Xdata), binaryc.predict(Xdata, baseline)



grid_fan = norm_var(0, grid_control[:, 0])

if USE_PUMP:

    grid_pump = norm_var(4, grid_control[:, 1])

binary_fan = norm_var(0, binary_control)



# Data with optimal control values

gridXdata, binaryXdata = np.copy(Xdata), np.copy(Xdata)

gridXdata[:, 0], binaryXdata[:, 0] = grid_control[:, 0], binary_control

if USE_PUMP:

    gridXdata[:, 4] = grid_control[:, 1]



grid_temp, binary_temp = est.predict(gridXdata), est.predict(binaryXdata)



plt.figure()

sgc = plt.plot(np.arange(len(grid_fan)), grid_fan, label='Grid Controller')

sbc = plt.plot(np.arange(len(grid_fan)), binary_fan, label='Binary Controller')

shc = plt.plot(np.arange(len(grid_fan)), Xfiltered[:, 0], label='Historic Control')

plt.ylabel('Fan Speed')

plt.xlabel('Time')

plt.sca(plt.gca().twinx())

sgt = plt.plot(np.arange(len(grid_fan)), norm_temp(grid_temp), ls=':', label='Optimal Temps')

sbt = plt.plot(np.arange(len(grid_fan)), norm_temp(binary_temp), ls=':', label='Binary Temps')

sht = plt.plot(np.arange(len(grid_fan)), norm_temp(Ydata), ls=':', label='Historic Temps')

sgt = plt.plot(np.arange(len(grid_fan)), norm_temp(grid_temp), ls=':', label='Optimal Temps', linewidth=2)

sbt = plt.plot(np.arange(len(grid_fan)), norm_temp(binary_temp), ls=':', label='Binary Temps', linewidth=2)

sht = plt.plot(np.arange(len(grid_fan)), norm_temp(Ydata), ls=':', label='Historic Temps', linewidth=2)

plt.ylabel('Temperature / K')

plots = sgc + sbc + shc + sgt + sbt + sht

labels = [p.get_label() for p in plots]

plt.legend(plots, labels)

plt.legend(plots, label)

ticks, labels = plt.xticks()

ticks = np.asarray([t for t in ticks if 0 <= t < len(tseries)])

plt.xticks(ticks, tseries[ticks.astype(int)], rotation=30)

plt.title('Control Performance ' + time.isoformat());

plt.show()

```



%% Cell type:code id: tags:



``` python

# Energy saving plots

XEnergydata = XEnergyfiltered

# New Xdata which is for plotting energy with the controller-set values

gridXEnergydata, binaryXEnergydata = np.copy(XEnergydata), np.copy(XEnergydata)

gridXEnergydata[:, 0], binaryXEnergydata[:, 0] = grid_temp, binary_temp

if USE_PUMP:

    gridXEnergydata[:, 7] = grid_control[:, 1]

# Model TempCondOut using chiller temperature model

chi_data = np.hstack((gridXEnergydata[:, 0:1], gridXEnergydata[:, 2:8]))

gridXEnergydata[:, 2] = est_chi.predict(chi_data)



grid_energy, binary_energy = est_energy.predict(gridXEnergydata), est_energy.predict(binaryXEnergydata)



plt.figure()

sgc = plt.scatter(np.arange(len(grid_energy)), norm_energy(grid_energy), label='Optimal Controller', marker='o')

sbc = plt.scatter(np.arange(len(grid_energy)), norm_energy(binary_energy), label='Binary Controller', marker='x')

sbc = plt.scatter(np.arange(len(grid_energy)), norm_energy(YEnergyfiltered), label='Historic', marker='+')

plt.xlabel('Time')

plt.ylabel('Power / W')

plt.title('Power usage over ' + time.isoformat())

plt.legend()

ticks, labels = plt.xticks()

ticks = np.asarray([t for t in ticks if 0 <= t < len(tseries)])

plt.xticks(ticks, tseries[ticks.astype(int)]);

plt.show()

```



%% Cell type:code id: tags:



``` python

# Comparison of control over a random sample

# where ambient temperature > 293 K = 20 C

temp_filter = Xtest[:, 1] > std_var(1, 273)

Xdata = Xtest[temp_filter]

Ydata = Ytest[temp_filter]



baseline = np.squeeze(std_temp(norm_var(3, Xdata[:, 3])))



grid_control, binary_control = gridc.predict(Xdata), binaryc.predict(Xdata, baseline)



grid_fan = norm_var(0, grid_control[:, 0])

if USE_PUMP:

    grid_pump = norm_var(4, grid_control[:, 1])

binary_fan = norm_var(0, binary_control)



# Data with values replaced by control set points

gridXdata, binaryXdata = np.copy(Xdata), np.copy(Xdata)

gridXdata[:, 0], binaryXdata[:, 0] = grid_control[:, 0], binary_control

if USE_PUMP:

    gridXdata[:, 4] = grid_control[:, 1]



grid_temp, binary_temp = est.predict(gridXdata), est.predict(binaryXdata)



plt.figure()

plt.subplot(121)

plt.hist(norm_var(0, grid_fan), label='Optimal Controller', histtype='step', density=False, bins=10, ls='-')

plt.hist(norm_var(0, binary_fan), label='Binary Controller', histtype='step', density=False, bins=10, ls=':')

plt.hist(norm_var(0, Xtest[:, 0]), label='Historic Control', histtype='step', density=False, bins=10, ls='-.')

plt.xlabel('PerFanFreq')

plt.ylabel('Frequency')

plt.ylim(0, 8000)

plt.legend()

plt.subplot(122)



plt.hist(norm_temp(grid_temp), label='Optimal Temps', histtype='step', ls='-')

plt.hist(norm_temp(binary_temp), label='Binary Temps', histtype='step', ls=':')

plt.hist(norm_temp(Ydata), label='Historic Temps', histtype='step', ls='-.')

plt.xlabel('Temperature / K')

plt.ylim(0, 8000)

plt.legend()

plt.suptitle('Control performance over test set')

plt.show();



print('Optimal\t\ttemp: mean {:.3f} std {:.3f}'.format(norm_temp(grid_temp.mean()), norm_temp(grid_temp).std()))

print('Binary\t\ttemp: mean {:.3f} std {:.3f}'.format(norm_temp(binary_temp.mean()), norm_temp(binary_temp).std()))

print('Historic\ttemp: mean {:.3f} std {:.3f}'.format(norm_temp(est.predict(Xdata).mean()), norm_temp(est.predict(Xdata)).std()))



print('Optimal\t\tfan: mean {:.3f} std {:.3f}'.format(norm_var(0, grid_fan.mean()), norm_var(0, grid_fan).std()))

print('Binary\t\tfan: mean {:.3f} std {:.3f}'.format(norm_var(0, binary_fan.mean()), norm_var(0, binary_fan).std()))

print('Historic\tfan: mean {:.3f} std {:.3f}'.format(norm_var(0, Xdata[:,0].mean()), norm_var(0, Xdata[:,0]).std()))

```



%% Cell type:code id: tags:



``` python

# Energy saving plots over a random sample

Xdata = XEnergytest[temp_filter]

Ydata = YEnergytest[temp_filter]



grid_temp, binary_temp = est.predict(gridXdata), est.predict(binaryXdata)



gridXEnergydata, binaryXEnergydata = np.copy(Xdata), np.copy(Xdata)

gridXEnergydata[:, 0], binaryXEnergydata[:, 0] = grid_temp, binary_temp

if USE_PUMP:

    gridXEnergydata[:, 7] = grid_control[:, 1]

# Model TempCondOut using chiller temperature model

chi_data = np.hstack((gridXEnergydata[:, 0:1], gridXEnergydata[:, 2:8]))

gridXEnergydata[:, 2] = est_chi.predict(chi_data)



grid_energy, binary_energy = est_energy.predict(gridXEnergydata), est_energy.predict(binaryXEnergydata)



plt.figure()

sgc = plt.hist(norm_energy(grid_energy), label='Grid Controller', histtype='step', ls='-')

sbc = plt.hist(norm_energy(binary_energy), label='Binary Controller', histtype='step', ls=':')

shc = plt.hist(norm_energy(Ydata), label='Historic', histtype='step', ls='-.')

sgc = plt.hist(norm_energy(grid_energy), label='Grid Controller', histtype='step', ls='-', lw=2)

sbc = plt.hist(norm_energy(binary_energy), label='Binary Controller', histtype='step', ls=':', lw=2)

shc = plt.hist(norm_energy(Ydata), label='Historic', histtype='step', ls='-.', lw=2)

plt.xlabel('Power / W')

plt.ylabel('Frequency')

plt.title('Power usage over test set')

plt.legend()

plt.show()



print('Optimal\t\tpower: mean {:.3f} std {:.3f}'.format(norm_energy(grid_energy).mean(), norm_energy(grid_energy).std()))

print('Binary\t\tpower: mean {:.3f} std {:.3f}'.format(norm_energy(binary_energy).mean(), norm_energy(binary_energy).std()))

print('Historic\tpower: mean {:.3f} std {:.3f}'.format(norm_energy(est_energy.predict(Xdata)).mean(), norm_energy(Ydata).std()))

```

src/Model.ipynb→src/Models.ipynb

+375 −329

File changed and moved.

Preview size limit exceeded, changes collapsed.

src/RL.ipynb

0 → 100644
+140 −0
Original line number Diff line number Diff line 
%% Cell type:code id: tags:



``` python

%matplotlib notebook

%reload_ext autoreload

%autoreload 2



import datetime

from os import path, environ

import pickle

import warnings



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, GridSearchCV, KFold

import torch

import torch.nn as nn



from utils import contiguous_sequences



# 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')

bin_path = './bin/'

```



%% Cell type:code id: tags:



``` python

# Data selection 'all' or 'chiller_on' or 'fan_on'

MODE = 'chiller_on'

# Read pre-processed data:

# Pytorch uses float32 as default type for weights etc,

# so input data points are also read in the same type.

df = pd.read_csv(chiller_file, index_col='Time',

                 parse_dates=['Time'], dtype=np.float32)

df.dropna(inplace=True)

if MODE == 'chiller_on':

    df = df[df['PowChi'] != 0.]

if MODE == 'fan_on':

    df = df[(df['PerFreqFanA'] != 0.) | df['PerFreqFanB'] != 0.]

print(len(df), 'Records')

```



%% Cell type:markdown id: tags:



## Environment model



%% Cell type:code id: tags:



``` python

inputs = ['PerFreqFans', 'PowConP', 'TempCondOut', 'TempAmbient', 'TempWetBulb']

outputs = ['PowChi', 'PowFans', 'TempCondOut']



df_in = pd.DataFrame(columns=inputs, index=df.index)

df_in['PerFreqFans'] = (df['PerFreqFanA'] + df['PerFreqFanB']) / 2

df_in[inputs[1:]] = df[inputs[1:]]



df_out = pd.DataFrame(columns=outputs, index=df.index)

df_out['PowChi'] = df['PowChi']

df_out['PowFans'] = df['PowFanA'] + df['PowFanB']

df_out['TempCondOut'] = df['TempCondOut']



idx_list = contiguous_sequences(df.index, pd.Timedelta(5, unit='min'))



dfs_in, dfs_out = [], []

for idx in idx_list:

    dfs_in.append(df_in.loc[idx[:-1]])

    dfs_out.append(df_out.loc[idx[1:]])



df_in = pd.concat(dfs_in, sort=False)

df_out = pd.concat(dfs_out, sort=False)



print('{:6d} time series'.format(len(dfs_in)))

print('{:6d} total rows'.format(len(df_in)))

```



%% Cell type:markdown id: tags:



### Fully-connected network



%% Cell type:code id: tags:



``` python

std_in, std_out = StandardScaler(), StandardScaler()

net = MLPRegressor(hidden_layer_sizes=(32, 32, 32),

                   activation='relu',

                   solver='adam',

                   verbose=True,

                   early_stopping=True)

est = Pipeline([('std', std_in), ('net', net)])



with warnings.catch_warnings():

    warnings.simplefilter('ignore', category=FutureWarning)

    est.fit(df_in, std_out.fit_transform(df_out))



```



%% Cell type:code id: tags:



``` python

save = {

    'loss': est['net'].loss_,

    'estimator': est,

    'output_norm': std_out,

    'inputs': inputs,

    'outputs': outputs

}

with open(path.join(bin_path, 'chiller_model_nn'), 'wb') as f:

    pickle.dump(save, f)

```



%% Cell type:code id: tags:



``` python

test_in, test_out = dfs_in[2], dfs_out[2]

pred = pd.DataFrame(std_out.inverse_transform(est.predict(test_in)),

                    index=test_out.index, columns=test_out.columns)



test_out.plot(grid=True)

pred.plot(grid=True)

```



%% Cell type:code id: tags:



``` python

# df_out.loc[idx_list[0]]

df_out[idx_list[0]]

```



%% Cell type:code id: tags:



``` python

df_out.head()

```

src/Relationships.ipynb

+1 −1

File changed.

Contains only whitespace changes.

src/controller.py

+0 −2
Original line number Diff line number Diff line
@@ -8,8 +8,6 @@ from sklearn.base import BaseEstimator 
import numpy as np

from scipy.optimize import fmin, minimize



from plotting import model_surface







class GridSearchController(BaseEstimator):