Commit 85dac0e6 authored by hazrmard's avatar hazrmard
Browse files

RL control with condenser water temp as setpoint; 

Renamed RL.ipynb > RL-Cooling Tower where fan speeds are setpoint,
RL-condenser.ipynb has simulations where condenser setpoint is action,
Changed model_surface to accept a plot function instead of a more specific BaseEstimator instance,
Fixes to docs (not yet complete)
parent 0927a5dc

docs/index.md

+4 −4
Changes for docs/index.md: 4 added lines, 4 removed lines.
Original line number Diff line number Diff line
@@ -26,16 +26,16 @@ The documentation is divided into a discussion of background concepts in physics 


4. [Chillers - Cooling towers](3-cooling-tower.md)



6. [Preprocessing](4-preprocessing.md)

5. Data



5. [Dataset and system description](5-dataset.md)

    a. [V1 Data](./datasets/v1/dataset.md)



    b. [V1 Preprocessing](./datasets/v1/preprocessing.md)



6. [Trends](6-trends.md)



7. [Relationships](7-relationships.md)



8. [Models](8-models.md)





## Installation

src/Models.ipynb

+4 −0
Changes for src/Models.ipynb: 4 added lines, 0 removed lines.
Original line number Diff line number Diff line 
%% Cell type:markdown id: tags:



Uses `v1` dataset.



%% 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

# Time at which to draw sample plots

# 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]

```



%% 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: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)],

    '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. Chiller 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

# Searching parameter grid for best hyperparameters

param_grid = {

    '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(XEnergytrain, YEnergytrain)



est = grid_search.best_estimator_

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



ENERGY_MODELS['SingleMLP'] = est

pickle.dump(est, open('./bin/PowChi_{}.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/PowChi_{}.pickle'.format(MODE), 'rb'))

ENERGY_MODELS['SingleMLP'] = est

```



%% Cell type:code id: tags:



``` python

# Plot historic predictions

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

# Searching parameter grid for best hyperparameters

param_grid = {

    'hidden_layer_sizes': [(4, 4, 4), (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(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

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

Xfiltered = X[select]

Yfiltered = Y[select]



# 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:markdown id: tags:



# 4. Cooling tower fan power model



%% Cell type:code id: tags:



``` python

# Data

fan_cols = ('PerFreqFanA', 'PerFreqFanB')

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

target_cols = ('PowFanA', 'PowFanB')

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

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

YPower = df.loc[:,target_cols].mean(axis=1)[:, None]

ScalerXPower, ScalerYPower = StandardScaler().fit(XPower), StandardScaler().fit(YPower)

XPower, YPower = ScalerXPower.transform(XPower), np.squeeze(ScalerYPower.transform(YPower))

XPowertrain, XPowertest, YPowertrain, YPowertest = train_test_split(XPower, YPower, test_size=0.1)



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

norm_power = lambda v: ScalerYPower.inverse_transform(v)

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



POWER_MODELS = {}

```



%% 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)],

    '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(XPowertrain, YPowertrain)



est = grid_search.best_estimator_

print('Test score:', est.score(XPowertest, YPowertest))



POWER_MODELS['SingleMLP'] = est

pickle.dump(est, open('./bin/power_{}.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/power_{}.pickle'.format(MODE), 'rb'))

POWER_MODELS['SingleMLP'] = est

```



%% Cell type:code id: tags:



``` python

# Plot historic predictions

Xfiltered = XPower[select]

Yfiltered = YPower[select]



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

            ScalerYPower.inverse_transform(Yfiltered),

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

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

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

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

plt.title('Fan 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, var_idx = 'PerFanFreq', 0



vary_range = (std_power_var(var_idx, 0),

              std_power_var(var_idx, 1))

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

                        (vary_range,), (20,))

ax = plot_surface(tseries,

                  norm_power_var(var_idx, y),

                  norm_power(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('Fan 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 = XPower[select], YPower[select]

var, var_idx = ('PerFreqFan', 'TempAmbient'), (0, 1)

vary_range = ((std_power_var(var_idx[0], 0), std_power_var(var_idx[0], 1)),

              (std_power_var(var_idx[1], 273), std_power_var(var_idx[1], 300)))



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

                        vary_range, (10, 10))

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

                  norm_power_var(var_idx[1], y),

                  norm_power(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('Fan Power Extrapolation ' + dt.isoformat());

```

src/RL-Condenser.ipynb

0 → 100644
+668 −0

File added.

Preview size limit exceeded, changes collapsed.

src/RL.ipynb→src/RL-Cooling Tower.ipynb

+5 −1
Changes for src/RL-Cooling Tower.ipynb: 5 added lines, 1 removed line.
Original line number Diff line number Diff line 
%% Cell type:markdown id: tags:



Uses `v1` dataset.



%% Cell type:code id: tags:



``` python

%matplotlib notebook

%reload_ext autoreload

%autoreload 2



import datetime

import sys

from os import path, environ

import pickle

import warnings



sys.path.insert(0, path.abspath('../../../pyTorchBridge/'))



import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

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 torch.optim as optim

from pytorchbridge import TorchEstimator

from tqdm.auto import tqdm, trange



from utils import contiguous_sequences

from preprocessing.thermo import k2f



# 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



### Cooling tower + chiller



%% Cell type:code id: tags:



``` python

inputs = ['PerFreqFans', 'PowConP', 'TempCondOut', 'TempAmbient', 'TempWetBulb', 'TempEvapIn', 'TempEvapOut', 'FlowEvap']

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

lag = 0    # 0, 1, 2, 3, 4, ...



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']

df_out['TempCondIn'] = df['TempCondIn']



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



# Create dataframes of contiguous sequences with a delay

# of 1 time unit to indicate causality input -> outputs

dfs_in, dfs_out = [], []

for idx in idx_list:

    dfs_in.append(df_in.loc[idx[:-lag if lag > 0 else None]])

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



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=(64, 32, 32, 16),

                   activation='tanh',

                   solver='adam',

                   verbose=True,

                   early_stopping=True,

                   learning_rate_init=1e-3)

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:markdown id: tags:



### Recurrent network



%% Cell type:code id: tags:



``` python

std_in, std_out = StandardScaler(), StandardScaler()

class MyModule(nn.Module):

    def __init__(self):

        super().__init__()

        self.lstm = nn.LSTM(input_size=5, hidden_size=3, num_layers=1, batch_first=True)

    def forward(self, x, h0=None):

        output, _ = self.lstm(x, h0)

        return output



module = MyModule()

loss = nn.MSELoss()

opt = optim.Adam(module.parameters(), lr=1e-2)

net = TorchEstimator(module, opt, loss, 30, verbose=True, batch_size=32)



std_out.fit(df_out)

std_in.fit(df_in)

X = [torch.as_tensor(std_in.transform(df.values)) for df in dfs_in]

Y = [torch.as_tensor(std_out.transform(df.values)) for df in dfs_out]

net.fit(nn.utils.rnn.pad_sequence(X, batch_first=True),

        nn.utils.rnn.pad_sequence(Y, batch_first=True))

```



%% Cell type:code id: tags:



``` python

# Save model

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:

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

    pickle.dump(save, f)

```



%% Cell type:code id: tags:



``` python

# load model

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

    save = pickle.load(f, fix_imports=False)

    est = save['estimator']

    std_out = save['output_norm']

```



%% Cell type:code id: tags:



``` python

# Visualize model predictions

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.loc[:, ('TempCondOut', 'TempCondIn')].plot(grid=True)

pred.loc[:, ('TempCondOut', 'TempCondIn')].plot(grid=True)

```



%% Cell type:markdown id: tags:



## RL Environment



%% Cell type:code id: tags:



``` python

# Make wrapper for cooling tower such that outputs are normalized

# i.e. in physical units instead of being 0 mean and 1 variance.

from cooling_tower import CoolingTower



externalvars = ('PowConP', 'TempAmbient', 'TempWetBulb', 'TempEvapIn', 'TempEvapOut', 'FlowEvap')

externalvals = [df.loc[:, externalvars].values for df in dfs_in]

is_recurrent = False

pump_control = False



class ESBModel:



    def __init__(self, estimator, transformer):

        self.estimator = estimator

        self.transformer = transformer



    def predict(self, x):

        return self.transformer.inverse_transform(self.estimator.predict(x))





esb = CoolingTower(ESBModel(est, std_out), is_recurrent, externalvals, pump_control)

```



%% Cell type:code id: tags:



``` python

# Visualize environment episode

done = False

states = []

power = []

esb.reset()

while not done:

    state, _, done, info = esb.step(esb.action_space.sample())

    states.append(state)

    power.append(info.get('powchi') + info.get('powfans'))

esb.reset()



states = np.asarray(states)

power = np.asarray(power)

plt.subplot(2,1,1)

plt.plot(states[:, 0], label='PowConP')

plt.plot(power, label='Total Power')

plt.legend()

plt.subplot(2,1,2)

plt.plot(states[:, 1], label='TempCondOut')

plt.plot(states[:, 2], label='TempAmbient')

plt.plot(states[:, 3], label='TempWetBulb')

plt.legend()

```



%% Cell type:markdown id: tags:



## RL Control



%% Cell type:code id: tags:



``` python

import tensorflow as tf

tf.compat.v1.logging.set_verbosity(tf.compat.v1.logging.ERROR)



from stable_baselines import PPO2

from stable_baselines.common.vec_env import DummyVecEnv

from stable_baselines.common.policies import MlpPolicy, LstmPolicy



class CT(CoolingTower):

    def reward(self, t, state: np.ndarray, action: np.ndarray, nstate: np.ndarray,

               locals: dict) -> float:

        powchi = (locals.get('powchi') - 5300) / 415700.

        powfans = locals.get('powfans') / 21230

        tempcondin = (locals.get('tempcondin') - 281.7) / 25.15

        tempcondout = (locals.get('tempcondout') - 290.) / 20.



        # return - (0.2 * powchi) - (0.8 * powfans)

        return - (0.9 * tempcondin) - (0.1 * powfans)

        # return - (0.1 * powchi) - (0.9 * tempcondin)



esb_vec = DummyVecEnv([lambda: CT(ESBModel(est, std_out), is_recurrent, externalvals, pump_control) \

                       for _ in range(10)])

agent = PPO2(MlpPolicy, esb_vec, verbose=1, learning_rate=25e-4)

agent.learn(150000, log_interval=10)

```



%% Cell type:code id: tags:



``` python

# seqidx = np.random.randint(len(dfs_in))

seqidx = 21 # August months

simulate_hist = False  # Whether to use raw output data, or simulate it through historical actions



# indexing histories after 1st element because simulated trajectories

# are recorded after initial state (> 0), so lengths are equal

act_hist = dfs_in[seqidx].loc[:, ('PerFreqFans', 'PowConP')].values[lag:]

ext = dfs_in[seqidx].loc[:, externalvars].values



# Get baseline by running historic actions through environment:

if simulate_hist:

    esb.reset(external=ext)

    done = False

    pow_hist_chi, pow_hist_fan, temp_hist = [], [], []

    t = 0

    while not done:

        action = act_hist[t, :2] if pump_control else act_hist[t, :1]

        _, _, done, info = esb.step(action)

        pow_hist_fan.append(info.get('powfans'))

        pow_hist_chi.append(info.get('powchi'))

        temp_hist.append(info.get('tempcondin'))

        t += 1

else:

    pow_hist_chi = dfs_out[seqidx]['PowChi'].values

    pow_hist_fans = dfs_out[seqidx]['PowFans'].values

    temp_hist = dfs_out[seqidx]['TempCondIn'].values

```



%% Cell type:code id: tags:



``` python

pfan, pchi, act, rewards, temp = [], [], [], [], []



# run multiple trials over same period for stochastic policy

for trial in range(10):

    state = esb.reset(external=ext)

    done = False

    pfan.append([])

    pchi.append([])

    act.append([])

    rewards.append([])

    temp.append([])

    while not done:

        action = agent.predict(state)[0]

        state, reward, done, info = esb.step(action)

        act[-1].append(action)

        pfan[-1].append(info.get('powfans'))

        pchi[-1].append(info.get('powchi'))

        rewards[-1].append(reward)

        temp[-1].append(info.get('tempcondin'))



# get std_dev and mean of metrics

std_pfan = np.std(pfan, axis=0, keepdims=False)

std_pchi = np.std(pchi, axis=0, keepdims=False)

std_act = np.std(act, axis=0, keepdims=False)

std_rewards = np.std(rewards, axis=0, keepdims=False)

std_temp = np.std(temp, axis=0, keepdims=False)



pfan = np.mean(pfan, axis=0, keepdims=False)

pchi = np.mean(pchi, axis=0, keepdims=False)

act = np.mean(act, axis=0, keepdims=False)

rewards = np.mean(rewards, axis=0, keepdims=False)

temp = np.mean(temp, axis=0, keepdims=False)

```



%% Cell type:code id: tags:



``` python

tempf = k2f(temp)  # convert to Farenheit

std_tempf = std_temp * 1.8

temp_histf = k2f(np.asarray(temp_hist))



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

plt.subplot(4,1,1)

plt.title('Fan Power (Average RL {:.0f}W vs Historical {:.0f}W)'\

          .format(np.mean(pfan), np.mean(pow_hist_fan)))

plt.plot(pfan, 'b:', label='RL.Fan')

plt.fill_between(np.arange(len(pfan)), pfan+std_pfan, pfan-std_pfan, color='b', alpha=0.3)

plt.plot(pow_hist_fan, 'r:', label='Historical.Fan')

plt.ylim(bottom=0)

plt.legend()



plt.subplot(4,1,2)

plt.title('Chiller Power (Average RL {:.0f}W vs Historical {:.0f}W)'\

          .format(np.mean(pchi), np.mean(pow_hist_chi)))

plt.plot(pchi, 'b:', label='RL.Chiller')

plt.fill_between(np.arange(len(pchi)), pchi+std_pchi, pchi-std_pchi, color='b', alpha=0.3)

plt.plot(pow_hist_chi, 'r:', label='Historical.Chiller')

plt.ylim(bottom=0)

plt.legend()



plt.subplot(4,1,3)

plt.title('Fan Control (Average RL {:.2f} vs Historical {:.2f})'\

          .format(np.mean(act[:, 0]), np.mean(act_hist[:, 0])))

plt.plot(act[:, 0], 'b:', label='RL.Fan')

plt.fill_between(np.arange(len(act[:, 0])), act[:, 0]+std_act[:, 0], act[:, 0]-std_act[:, 0], color='b', alpha=0.3)

plt.plot(act_hist[:, 0], 'r:', label='Historical.Fan')

plt.ylim(top=1.05)

plt.legend()





plt.subplot(4,1,4)

plt.title('Output Temperature (Average RL {:.1f}F vs Historical {:.1f}F)'\

          .format(np.mean(tempf), np.mean(temp_histf)))

plt.plot(tempf, 'b:', label='RL.Temp')

plt.fill_between(np.arange(len(tempf)), tempf+std_tempf, tempf-std_tempf, color='b', alpha=0.3)

plt.plot(temp_histf, 'r:', label='Historical.Temp')

plt.legend()



plt.tight_layout()

```



%% Cell type:code id: tags:



``` python

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

plt.plot(k2f(dfs_in[seqidx]['TempAmbient'].values), label='Ambient Temp')

plt.plot(k2f(dfs_in[seqidx]['TempWetBulb'].values), label='WetBulb Temp')

plt.legend()

plt.ylabel('Temperature /F')

plt.title('Environmental Conditions')

plt.tight_layout()

plt.show()

```

src/Trends.ipynb

+4 −1
Changes for src/Trends.ipynb: 4 added lines, 1 removed line.
Original line number Diff line number Diff line 
%% Cell type:markdown id: tags:



Uses `v1` dataset.



%% Cell type:code id: tags:



``` python

%matplotlib notebook

%reload_ext autoreload

%autoreload 2



from datetime import datetime

from os import path, environ



import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from matplotlib.animation import FuncAnimation

import ipyvolume as ipv



from preprocess import POW_FIELDS, TEMP_FIELDS

from thermo import CONSTANTS

from plotting import animate_dataframes

# 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:markdown id: tags:



# Trends



## Dataset



The dataset describes input/output powers, temperatures, and cooling by various components of a chiller with a water-cooled condenser and a cooling tower. More details are in `/docs/`.



After pre-processing, all measurements are in SI units (Kelvins, watts) except for the field `KWPerTon` where units are explicit in the measurement.



%% 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)

print('{} rows, {} columns'.format(len(df), len(df.columns)))

df.head()

```



%% Cell type:markdown id: tags:



## Efficiency metrics



### Coefficient of performance



$$

\begin{align*}

COP_{max cooling} &= \frac{T_{cold}}{T_{hot} - T_{cold}}  \\

COP &= \frac{\texttt{Energy Extracted}}{\texttt{Energy Input}}

\end{align*}

$$



%% Cell type:code id: tags:



``` python

# Measure maximum Coefficient of Performance (COP) for cooling

# for the ENTIRE chiller plant.

# Currently in the dataset, water flow rate for the condenser loop

# is not available - so COP for condenser/cooling tower cannot be

# calculated.

df['COPMax'] = df['TempEvapIn'] / (df['TempEvapIn'] - df['TempEvapOut'])



# Measured achieved COP

df['COP'] = df['PowCool'] / df['PowIn']

```



%% Cell type:code id: tags:



``` python

# Aggregating measurements across days by time of day

df = df[~df.index.duplicated(keep='first')]          # temporary fix to ignore duplicate timestamps

downsampled = df.asfreq('15T')                       # downsampling to 15 min for cleaner plot



# Grouping the data by week and taking means over each time

weekly = downsampled.resample('W') # iterator over (date week ends, dataframe for week)

wlabels = [ts.strftime('Week ending: %Y-%b-%d') for ts in weekly.groups.keys()]

means_by_week = []                 # list of dataframes of timestamp means over each week

for date, subframe in weekly:

    grouped = subframe.groupby(subframe.index.time)

    week_mean = grouped.mean()

    week_mean.set_index(pd.to_datetime(week_mean.index, format='%H:%M:%S').time, inplace=True)

    means_by_week.append(week_mean)

```



%% Cell type:markdown id: tags:



## Plots



### COP



%% Cell type:code id: tags:



``` python

series = ['COP', 'COPMax']

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

f.suptitle('COP (Entire plant) - Daily Averages')

ax = f.add_subplot(111)

anim = animate_dataframes(frames=means_by_week, ax=ax, lseries=series, labels=wlabels,

                          ylim=(1e-2, 500), xlim=(min(df.index.time), max(df.index.time)),

                          xlabel='Time', ylabel='COP',

                          yscale='log', anim_args={'repeat':False, 'blit':True})

anim.save(path.join(plot_path, '6-COP.mp4'))

```



%% Cell type:markdown id: tags:



### Temperatures



%% Cell type:code id: tags:



``` python

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

f.suptitle('Temperature - Daily Averages')

ax = f.add_subplot(111)

anim = animate_dataframes(frames=means_by_week, ax=ax, lseries=TEMP_FIELDS, labels=wlabels,

                          ylim=(273, 313), xlim=(min(df.index.time), max(df.index.time)),

                          xlabel='Time', ylabel='Temperature (K)',

                          anim_args={'repeat':False, 'blit':True})

anim.save(path.join(plot_path, '6-temps.mp4'))

```



%% Cell type:markdown id: tags:



### Power vs. Cooling



Input power vs. total cooling done by the evaporator.



%% Cell type:code id: tags:



``` python

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

f.suptitle('Input power vs. Cooling')

series = ['PowIn', 'PowCool']

ax = f.add_subplot(111)

anim = animate_dataframes(frames=means_by_week, ax=ax, lseries=series, labels=wlabels,

                          ylim=None, xlim=(min(df.index.time), max(df.index.time)),

                          xlabel='Time', ylabel='Power (watts)',

                          anim_args={'repeat':False, 'blit':True})

anim.save(path.join(plot_path, '6-pow-vs-cooling.mp4'))

```



%% Cell type:markdown id: tags:



### Condenser/cooling tower cycle power



%% Cell type:code id: tags:



``` python

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

f.suptitle('Cooling tower power vs cooling')

ax = f.add_subplot(111)

lseries = ['PowFanA', 'PowFanB', 'PowConP']

rseries = ['TempCondOut', 'TempCondIn', 'TempAmbient', 'TempWetBulb']

anim = animate_dataframes(frames=means_by_week, ax=ax, lseries=lseries, rseries=rseries,

                          labels=wlabels,

                          ylim=((0,np.nanmax(df[lseries].values)), (273, 313)),

                          xlim=(min(df.index.time), max(df.index.time)),

                          xlabel='Time', ylabel=('Power (watts)', 'Temp(K)'),

                          anim_args={'repeat':False, 'blit':True})

anim.save(path.join(plot_path, '6-cooling-tower-pow.mp4'))

```



%% Cell type:markdown id: tags:



### Fan speed (%) vs fan power



%% Cell type:code id: tags:



``` python

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

f.suptitle('Fan power (w) and fan speed (%)')

ax = f.add_subplot(111)

lseries = ['PowFanA', 'PowFanB']

rseries = ['PerFreqFanA', 'PerFreqFanB']

anim = animate_dataframes(frames=means_by_week, ax=ax, lseries=lseries, rseries=rseries,

                          labels=wlabels,

                          ylim=((0,np.nanmax(df[lseries].values)), (0, 1)),

                          xlim=(min(df.index.time), max(df.index.time)),

                          xlabel='Time', ylabel=('Power (Watts)', 'Percentage speed (%)'),

                          anim_args={'repeat':False, 'blit':True})

anim.save(path.join(plot_path, '6-fan-power-speed.mp4'))

```



%% Cell type:markdown id: tags:



### Distribution of fan power signals



%% Cell type:code id: tags:



``` python

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

ax = f.add_subplot(111)

f.suptitle('Distribution of Fan power signals')

df['PerFreqFanA'].plot.hist(bins=20, legend=True, ax=ax, histtype='step')

df['PerFreqFanB'].plot.hist(bins=20, legend=True, ax=ax, histtype='step')

plt.savefig(path.join(plot_path, '6-fan-power-hist.png'))

```