Commit 52a1f7d4 authored by Ibrahim Ahmed's avatar Ibrahim Ahmed
Browse files

baseline_control/controller: Added randomness to actions

parent cee659a5
Loading
Loading
Loading
Loading
+1 −0
Changes for .gitignore: 1 added line, 0 removed lines.
Original line number Diff line number Diff line
credentials.txt
questions.txt
src/Scratch.ipynb
*.xlsx
*.pptx
*.pdf
+49 −4
Changes for src/Baseline-Condenser.ipynb: 49 added lines, 4 removed lines.
Original line number Diff line number Diff line
%% Cell type:code id: tags:

``` python
%matplotlib inline
%reload_ext autoreload
%autoreload 2

import datetime
import sys
from os import path, environ
import pickle
import warnings

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 tqdm.auto import tqdm, trange

from utils import contiguous_sequences
from plotting import model_surface, plot_surface
from systems import Condenser
from baseline_control import SimpleFeedbackController, FeedbackController

chiller_file_1 = path.join(environ['DATADIR'],
                         'EngineeringScienceBuilding',
                         '2422_ESB_HVAC_1.csv')
chiller_file_2 = path.join(environ['DATADIR'],
                         'EngineeringScienceBuilding',
                         '2841_ESB_HVAC_2.csv')

plot_path = path.join('..', 'docs', 'img')
bin_path = './bin/'
```

%% Cell type:markdown id: tags:

# Controller script functions
## Controller script demo

%% Cell type:code id: tags:

``` python
from datetime import datetime, timedelta
import pytz
from controller import make_arguments, get_settings, get_controller, update_controller, get_current_state, put_control_action
```

%% Cell type:code id: tags:

``` python
parser = make_arguments()
args = parser.parse_args(['-s', './local.ini'])
settings = get_settings(args)
settings['target'] = 'temperature'
settings
```

%% Cell type:code id: tags:

``` python
ctrl = get_controller(**settings)
update_controller(ctrl, **settings)
```

%% Cell type:code id: tags:

``` python
end = datetime.now(pytz.utc)
start = end - timedelta(minutes=10)
s = get_current_state(start, end, **settings)
s
```

%% Cell type:code id: tags:

``` python
s['TempCondIn'] = 68
ctrl = get_controller(**settings)
actions = []
feedbacks = []
temps = []
s['TempCondIn'] = 62.
s['TempWetBulb'] = 40.
for i in range(40):
    action, = ctrl.predict(s)
    feedbacks.append(ctrl.feedback(s))
    actions.append(action)
    temps.append(s['TempCondIn'])
    if i < 10:
        s['TempCondIn'] -= 1.
    elif i < 20:
        s['TempCondIn'] += 1.
    elif i < 30:
        if actions[-1] > actions[-2]:
            s['TempCondIn'] -= 1.
        else:
            s['TempCondIn'] += 1.
    elif i < 40:
        if actions[-1] > actions[-2]:
            s['TempCondIn'] += 1.
        else:
            s['TempCondIn'] -= 1.
```

%% Cell type:code id: tags:

``` python
action, = ctrl.predict(s)
action
plt.figure(figsize=(12,8))
# plt.imshow(np.zeros((1, 20)), aspect='auto', alpha=0.3)
plt.plot(temps, label='Temp /F')
plt.plot(actions, label='Setpoint /F')
plt.grid(which='both')
for line in (10,20,30):
    plt.axvline(x=line, color='black', ls=':')
plt.axhline(y=s['TempWetBulb'], label='WetBulb /F', color='red', ls='--')
plt.axhline(y=55, label='Action lower bound', color='blue', ls='--')
plt.legend(loc='upper left')

plt.text(2, 47, 'Increasing\n(Unresponsive)')
plt.text(12, 47, 'Decreasing\n(Unresponsive)')
plt.text(22, 47, 'Same direction')
plt.text(32, 47, 'Opposite direction')
plt.title('Controller response to different feedback behaviors')
plt.xlabel('Time')
plt.ylabel('Temperature /F')

plt.twinx()
plt.plot(feedbacks, 'g:', lw=3, label='feedback')
plt.legend(loc='upper right')
plt.ylabel('Feecback /F')
```

%% Cell type:markdown id: tags:

## Environment model

State variables (12):

`'TempCondIn', 'TempCondOut', 'TempEvapOut', 'PowChi', 'PowFanA', 'PowFanB', 'PowConP', 'TempEvapIn', 'TempAmbient', 'TempWetBulb', 'PressDiffEvap', 'PressDiffCond'`

Action variables (1):

`'TempCondInSetpoint'`

Output variables (3):

`'TempCondIn', 'TempCondOut', 'TempEvapOut', 'PowChi', 'PowFanA', 'PowFanB', 'PowConP'`

Model:

`[Action, State] --> [Output]`

%% Cell type:code id: tags:

``` python
# Choosing which chiller to use
chiller_file = chiller_file_2
# 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)
print('Original length: {} Records'.format(len(df)))
# # These fields were not populated until 2020-07-01, so leaving then out of analysis
# df.drop(['PowFanA', 'PowFanB', 'FlowCond', 'PowChiP', 'PerFreqConP', 'PowConP'], axis='columns', inplace=True)
df.drop(['FlowCond', 'PowChiP', 'PerFreqConP'], axis='columns', inplace=True)

df.dropna(inplace=True)
if MODE == 'chiller_on':
    df = df[df['RunChi'] != 0]
if MODE == 'fan_on':
    df = df[(df['RunFanA'] != 0.) | df['RunFanB'] != 0.]
print('Processed length: {} Records'.format(len(df)))
```

%% Cell type:code id: tags:

``` python
# load model
with open(path.join(bin_path, 'v2_condenser'), 'rb') as f:
    save = pickle.load(f, fix_imports=False)
    est_cond = save['estimator']
    std_out_cond = save['output_norm']
    statevars = save['statevars']
    actionvars = save['actionvars']
    inputs = save['inputs']
    outputs = save['outputs']
    lag = save['lag']
```

%% Cell type:markdown id: tags:

### Condenser Data

%% Cell type:code id: tags:

``` python
df_in = pd.DataFrame(columns=inputs, index=df.index)
df_in['TempCondInSetpoint'] = np.clip(df['TempWetBulb'] - 4, a_min=65, a_max=None)  # approach controller
df_in[inputs[1:]] = df[inputs[1:]]

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

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[:-max(lag) if max(lag) > 0 else None]])
    cols = []
    for l, c in zip(lag, outputs):
        window = slice(l, None if l==max(lag) else -(max(lag)-l))
        series = df_out[c].loc[idx[window]]
        cols.append(series.values)
        if l == min(lag): index = series.index
    dfs_out.append(pd.DataFrame(np.asarray(cols).T, index=index, columns=outputs))

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:

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

externalvars = ('TempEvapIn', 'TempAmbient', 'TempWetBulb', 'PressDiffEvap', 'PressDiffCond')
externalvals = [df.loc[:, externalvars] for df in dfs_in]

class InvTransformer:

    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 = Condenser(InvTransformer(est_cond, std_out_cond), externalvals)
```

%% 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'))
esb.reset()

states = np.asarray(states)
power = np.asarray(power)
plt.subplot(2,1,1)
plt.plot(power, label='Total Power')
plt.legend()
plt.subplot(2,1,2)
plt.plot(states[:, 0], label='TempCondIn')
plt.plot(states[:, 1], label='TempCondOut')
plt.plot(states[:, 2], label='TempEvapOut')
plt.plot(states[:, 4], label='TempEvapIn')
plt.plot(states[:, 5], label='TempAmbient')
plt.plot(states[:, 6], label='TempWetBulb')
plt.legend()
plt.show()
```

%% Cell type:markdown id: tags:

## Simple Feedback Control

%% Cell type:code id: tags:

``` python
longest_seq_idx = max(range(len(dfs_in)), key= lambda i: len(dfs_in[i]))
```

%% Cell type:code id: tags:

``` python
dfs_in[longest_seq_idx]
```

%% Cell type:code id: tags:

``` python
# seqidx = np.random.randint(len(dfs_in))
seqidx = longest_seq_idx
simulate_hist = True  # 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[:, ('TempCondInSetpoint')].values[1:, None]
ext = dfs_in[seqidx].loc[:, externalvars]

# Get baseline by running historic actions through environment:
if simulate_hist:
    esb.reset(external=ext, state0=dfs_in[seqidx].iloc[0, 1:].values)
    done = False
    pow_hist_chi, pow_hist_fan, temp_hist = [], [], []
    t = 0
    while not done:
        action = 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
# define agent
class Controller1(SimpleFeedbackController):

    def feedback(self, X):
        return -sum(X[3:7])  # PowChi, PowFanA, PowFanB, PowConP
        # return -X[0]

    def starting_action(self, X):
        return np.asarray([X[9] + 4]) # TempWetBulb

    def clip_action(self, u, X):
        u = super().clip_action(u, X)
        return np.clip(u, a_min=X[9], a_max=None)

class Controller2(FeedbackController):

    def feedback(self, X):
        return -X[3]  # PowChi

    def starting_action(self, X):
        return None
        # return X[9] + 4 # TempWetBulb



agent_fn = lambda: Controller1(bounds=((60., 80.),), stepsize=1)
# agent_fn = lambda: Controller2(bounds=((55., 90.),), kp=1., ki=0.2, kd=0.)
```

%% Cell type:code id: tags:

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

# run multiple trials over same period for stochastic policy
for trial in trange(1, leave=False):
    state = esb.reset(external=ext, state0=dfs_in[seqidx].iloc[0, 1:].values)
    agent = agent_fn()
    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
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(3,1,1)
plt.title('Chiller Power (Average {:.0f}kW vs Historical {:.0f}kW)'\
          .format(np.mean(pchi), np.mean(pow_hist_chi)))
plt.plot(pchi, 'b:', label='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(3,1,2)
plt.title('Setpoint Control (Average {:.2f} vs Historical {:.2f})'\
          .format(np.mean(act[:, 0]), np.mean(act_hist[:, 0])))
plt.plot(ext['TempAmbient'].values, 'g.', label='TempAmbient')
plt.plot(ext['TempWetBulb'].values, 'c.', label='TempWetBulb')
plt.plot(act[:, 0], 'b:', label='Setpoint')
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.Setpoint')
# plt.ylim(top=1.05)
plt.legend()


plt.subplot(3,1,3)
plt.title('Output Temperature (Average {:.1f}F vs Historical {:.1f}F)'\
          .format(np.mean(temp), np.mean(temp_hist)))
plt.plot(temp, 'b:', label='Temp')
plt.fill_between(np.arange(len(temp)), temp+std_temp, temp-std_temp, color='b', alpha=0.3)
plt.plot(temp_hist, 'r:', label='Historical.Temp')
plt.legend()

plt.tight_layout()
```

%% Cell type:code id: tags:

``` python
plt.figure(figsize=(8,3))
plt.plot(dfs_in[seqidx]['TempAmbient'].values, label='Ambient Temp')
plt.plot(dfs_in[seqidx]['TempWetBulb'].values, label='WetBulb Temp')
plt.legend()
plt.ylabel('Temperature /F')
plt.title('Environmental Conditions')
plt.tight_layout()
plt.show()
```
+56 −135
Changes for src/RL-Condenser.ipynb: 56 added lines, 135 removed lines.
Original line number Diff line number Diff line
%% 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 tqdm.auto import tqdm, trange

from utils import contiguous_sequences
from plotting import model_surface, plot_surface
from condenser import Condenser

chiller_file = path.join(environ['DATADIR'],
                         'EngineeringScienceBuilding',
                         '2422_ESB_HVAC.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.drop(['PowFanA', 'PowFanB', 'FlowCond', 'PowChiP', 'PerFreqConP', 'PowConP'], axis='columns', inplace=True)
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

State variables (8):

`'TempCondIn', 'TempCondOut', 'TempAmbient', 'TempWetBulb', 'TempEvapIn', 'TempEvapOut', 'PressDiffEvap', 'PressDiffCond'`

Action variables (1):

`'TempCondInSetpoint'`

%% Cell type:markdown id: tags:

### Condenser

%% Cell type:code id: tags:

``` python
envvars = ['TempCondIn', 'TempCondOut', 'TempAmbient', 'TempWetBulb', 'TempEvapIn', 'TempEvapOut', 'PressDiffEvap', 'PressDiffCond']
actionvars = ['TempCondInSetpoint']
inputs =  actionvars + envvars
outputs = ['PowChi', 'TempCondOut', 'TempCondIn']
lag = (1, 1, 1)    # 0, 1, 2, 3, 4, ...

df_in = pd.DataFrame(columns=inputs, index=df.index)
df_in['TempCondInSetpoint'] = np.clip(df['TempWetBulb'] - 4, a_min=65, a_max=None)  # approach controller
df_in[inputs[1:]] = df[inputs[1:]]

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

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[:-max(lag) if max(lag) > 0 else None]])
    cols = []
    for l, c in zip(lag, outputs):
        window = slice(l, None if l==max(lag) else -(max(lag)-l))
        series = df_out[c].loc[idx[window]]
        cols.append(series.values)
        if l == min(lag): index = series.index
    dfs_out.append(pd.DataFrame(np.asarray(cols).T, index=index, columns=outputs))

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

``` python
std_in_cond, std_out_cond = StandardScaler(), StandardScaler()
net = MLPRegressor(hidden_layer_sizes=(64, 32, 16),
                   activation='tanh',
                   solver='adam',
                   verbose=True,
                   early_stopping=True,
                   learning_rate_init=1e-3)
est_cond = Pipeline([('std', std_in_cond), ('net', net)])

with warnings.catch_warnings():
    warnings.simplefilter('ignore', category=FutureWarning)
    est_cond.fit(df_in, std_out_cond.fit_transform(df_out))
```

%% Cell type:code id: tags:

``` python
# Save model
save = {
    'loss': est_cond['net'].loss_,
    'estimator': est_cond,
    'output_norm': std_out_cond,
    'inputs': inputs,
    'outputs': outputs
}
with open(path.join(bin_path, 'env_condenser_condenser_nn'), 'wb') as f:
    pickle.dump(save, f)
```

%% Cell type:code id: tags:

``` python
# load model
with open(path.join(bin_path, 'env_condenser_condenser_nn'), 'rb') as f:
    save = pickle.load(f, fix_imports=False)
    est_cond = 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_cond.predict(test_in)),
                    index=test_out.index, columns=test_out.columns)

plt.subplot(2, 1, 1)
test_in.loc[:, ('TempCondIn')].plot(grid=True, style=':', label='TempCondIn')
test_out.loc[:, ('TempCondIn')].plot(grid=True, style=':', label='TempCondIn-Next')
pred.loc[:, ('TempCondIn')].plot(grid=True, style=':', label='TempCondIn-Next-Pred')
test_in.loc[:, ('TempCondOut')].plot(grid=True, style=':', label='TempCondOut')
test_out.loc[:, ('TempCondOut')].plot(grid=True, style=':', label='TempCondOut-Next')
pred.loc[:, ('TempCondOut')].plot(grid=True, style=':', label='TempCondOut-Next-Pred')
plt.legend()
plt.subplot(2, 1, 2)
ax4 = test_out.loc[:, ('PowChi')].plot(grid=True, label='PowChi')
ax5 = pred.loc[:, ('PowChi')].plot(grid=True, label='PowChi-Pred')
plt.legend()
```

%% Cell type:code id: tags:

``` python
point = df_in.loc['2019-07-01T1200-6'].values.reshape(1, -1)
x, y, z = model_surface(lambda x: std_out_cond.inverse_transform(est_cond.predict(x))[:,2],
                        X=point, vary_idx=(0, 3),vary_range=((65, 85), (75, 95)), vary_num=(20, 20))
ax = plot_surface(x,y,z, cmap=plt.cm.coolwarm)
ax.set_xlabel('TempCondInSetpoint')
ax.set_ylabel('TempAmbient')
ax.set_zlabel('TempCondIn-Next Cycle')
```

%% Cell type:code id: tags:

``` python
print(df_in.loc['2019-07-01T1200-6'])
```

%% Cell type:markdown id: tags:

### Cooling Tower

Note: not being used in the environment. The condenser model is now predicting the next timestep's `TempCondIn` as well.

%% Cell type:code id: tags:

``` python
envvars = ['TempCondOut', 'TempAmbient', 'TempWetBulb', 'TempEvapIn', 'TempEvapOut', 'PressDiffEvap', 'PressDiffCond']
actionvars = ['TempCondInSetpoint']
inputs =  actionvars + envvars
outputs = ['TempCondIn']
lag = 0    # 0, 1, 2, 3, 4, ...

df_in = pd.DataFrame(columns=inputs, index=df.index)
df_in['TempCondInSetpoint'] = df['TempWetBulb'] + 4  # approach controller
df_in[inputs[1:]] = df[inputs[1:]]

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

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

``` python
std_in_tower, std_out_tower = StandardScaler(), StandardScaler()
net = MLPRegressor(hidden_layer_sizes=(64, 32, 16),
                   activation='tanh',
                   solver='adam',
                   verbose=True,
                   early_stopping=True,
                   learning_rate_init=1e-3)
est_tower = Pipeline([('std', std_in_tower), ('net', net)])

with warnings.catch_warnings():
    warnings.simplefilter('ignore', category=FutureWarning)
    est_tower.fit(df_in, std_out_tower.fit_transform(df_out))
```

%% Cell type:code id: tags:

``` python
# Save model
save = {
    'loss': est_tower['net'].loss_,
    'estimator': est_tower,
    'output_norm': std_out_tower,
    'inputs': inputs,
    'outputs': outputs
}
with open(path.join(bin_path, 'env_condenser_tower_nn'), 'wb') as f:
    pickle.dump(save, f)
```

%% Cell type:code id: tags:

``` python
# load model
with open(path.join(bin_path, 'env_condenser_tower_nn'), 'rb') as f:
    save = pickle.load(f, fix_imports=False)
    est_tower = save['estimator']
    std_out_tower = 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_tower.inverse_transform(est_tower.predict(test_in)),
                    index=test_out.index, columns=test_out.columns)
test_in.loc[:, ('TempCondInSetpoint', 'TempCondOut')].plot(grid=True)
test_out.loc[:, ('TempCondIn')].plot(grid=True, style=':')
pred.loc[:, ('TempCondIn')].plot(grid=True, style=':')
plt.legend()
```

%% 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 = ('TempAmbient', 'TempWetBulb', 'TempEvapIn', 'TempEvapOut', 'PressDiffEvap', 'PressDiffCond')
externalvals = [df.loc[:, externalvars] for df in dfs_in]

class InvTransformer:

    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 = Condenser(InvTransformer(est_cond, std_out_cond), externalvals)
```

%% 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'))
esb.reset()

states = np.asarray(states)
power = np.asarray(power)
plt.subplot(2,1,1)
plt.plot(power, label='Total Power')
plt.legend()
plt.subplot(2,1,2)
plt.plot(states[:, 0], label='TempCondIn')
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)
        return -powchi

esb_vec = DummyVecEnv([lambda: Condenser(InvTransformer(est_cond, std_out_cond), externalvals) \
                       for _ in range(4)])
agent = PPO2(MlpPolicy, esb_vec, verbose=1, learning_rate=1e-3)
agent.learn(50000, log_interval=10)
```

%% Output

    -------------------------------------
    | approxkl           | 0.0043457365 |
    | clipfrac           | 0.04736328   |
    | explained_variance | 0.807        |
    | fps                | 490          |
    | n_updates          | 1            |
    | policy_entropy     | 1.4178611    |
    | policy_loss        | -0.004554482 |
    | serial_timesteps   | 128          |
    | time_elapsed       | 0            |
    | total_timesteps    | 512          |
    | value_loss         | 0.003511376  |
    -------------------------------------
    ---------------------------------------
    | approxkl           | 0.0068328762   |
    | clipfrac           | 0.100097656    |
    | explained_variance | 0.691          |
    | fps                | 1517           |
    | n_updates          | 10             |
    | policy_entropy     | 1.3978771      |
    | policy_loss        | -0.0037712269  |
    | serial_timesteps   | 1280           |
    | time_elapsed       | 3.76           |
    | total_timesteps    | 5120           |
    | value_loss         | 0.000100123114 |
    ---------------------------------------
    --------------------------------------
    | approxkl           | 0.00047332034 |
    | clipfrac           | 0.0           |
    | explained_variance | 0.832         |
    | fps                | 1495          |
    | n_updates          | 20            |
    | policy_entropy     | 1.4107119     |
    | policy_loss        | -0.0008063864 |
    | serial_timesteps   | 2560          |
    | time_elapsed       | 7.22          |
    | total_timesteps    | 10240         |
    | value_loss         | 9.3943345e-06 |
    --------------------------------------
    ---------------------------------------
    | approxkl           | 0.0032120398   |
    | clipfrac           | 0.041015625    |
    | explained_variance | -4.33          |
    | fps                | 1513           |
    | n_updates          | 30             |
    | policy_entropy     | 1.4007572      |
    | policy_loss        | -0.00067680766 |
    | serial_timesteps   | 3840           |
    | time_elapsed       | 10.6           |
    | total_timesteps    | 15360          |
    | value_loss         | 6.4343025e-05  |
    ---------------------------------------
    --------------------------------------
    | approxkl           | 0.005868512   |
    | clipfrac           | 0.08886719    |
    | explained_variance | 0.3           |
    | fps                | 1513          |
    | n_updates          | 40            |
    | policy_entropy     | 1.3987796     |
    | policy_loss        | -0.0045702755 |
    | serial_timesteps   | 5120          |
    | time_elapsed       | 14            |
    | total_timesteps    | 20480         |
    | value_loss         | 7.2725215e-06 |
    --------------------------------------
    --------------------------------------
    | approxkl           | 0.0020880182  |
    | clipfrac           | 0.022949219   |
    | explained_variance | 0.426         |
    | fps                | 1517          |
    | n_updates          | 50            |
    | policy_entropy     | 1.3999249     |
    | policy_loss        | -0.0010632982 |
    | serial_timesteps   | 6400          |
    | time_elapsed       | 17.4          |
    | total_timesteps    | 25600         |
    | value_loss         | 3.131517e-05  |
    --------------------------------------
    --------------------------------------
    | approxkl           | 0.0024089722  |
    | clipfrac           | 0.028320312   |
    | explained_variance | 0.362         |
    | fps                | 1531          |
    | n_updates          | 60            |
    | policy_entropy     | 1.3741724     |
    | policy_loss        | -0.0027649545 |
    | serial_timesteps   | 7680          |
    | time_elapsed       | 20.8          |
    | total_timesteps    | 30720         |
    | value_loss         | 7.0386354e-06 |
    --------------------------------------
    --------------------------------------
    | approxkl           | 0.006138832   |
    | clipfrac           | 0.08935547    |
    | explained_variance | 0.81          |
    | fps                | 1522          |
    | n_updates          | 70            |
    | policy_entropy     | 1.3747615     |
    | policy_loss        | -0.004558214  |
    | serial_timesteps   | 8960          |
    | time_elapsed       | 24.2          |
    | total_timesteps    | 35840         |
    | value_loss         | 1.7150476e-05 |
    --------------------------------------
    --------------------------------------
    | approxkl           | 0.00460754    |
    | clipfrac           | 0.061035156   |
    | explained_variance | 0.815         |
    | fps                | 1513          |
    | n_updates          | 80            |
    | policy_entropy     | 1.3724912     |
    | policy_loss        | -0.0011459764 |
    | serial_timesteps   | 10240         |
    | time_elapsed       | 27.6          |
    | total_timesteps    | 40960         |
    | value_loss         | 1.1484011e-05 |
    --------------------------------------
    ---------------------------------------
    | approxkl           | 0.0054067136   |
    | clipfrac           | 0.067871094    |
    | explained_variance | 0.165          |
    | fps                | 1526           |
    | n_updates          | 90             |
    | policy_entropy     | 1.3891732      |
    | policy_loss        | -0.00041625166 |
    | serial_timesteps   | 11520          |
    | time_elapsed       | 31             |
    | total_timesteps    | 46080          |
    | value_loss         | 2.0953019e-05  |
    ---------------------------------------

    <stable_baselines.ppo2.ppo2.PPO2 at 0x19a14de7a88>

%% Cell type:code id: tags:

``` python
# seqidx = np.random.randint(len(dfs_in))
seqidx = 21 # May
simulate_hist = True  # 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[:, ('TempCondInSetpoint')].values[1:, None]
ext = dfs_in[seqidx].loc[:, externalvars]

# 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, :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
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('Setpoint Control (Average RL {:.2f} vs Historical {:.2f})'\
          .format(np.mean(act[:, 0]), np.mean(act_hist[:, 0])))
plt.plot(act[:, 0], 'b:', label='RL.Setpoint')
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.Setpoint')
# 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(temp), np.mean(temp_hist)))
plt.plot(temp, 'b:', label='RL.Temp')
plt.fill_between(np.arange(len(temp)), temp+std_temp, temp-std_temp, color='b', alpha=0.3)
plt.plot(temp_hist, '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()
```
+27 −12
Changes for src/baseline_control.py: 27 added lines, 12 removed lines.
Original line number Diff line number Diff line
"""
Defines controller classes implementing various approaches.
Defines controller classes implementing various approaches. The controllers implemented
here do not use machine learning to make decisions. Although they may optionally employ
machine-learned models for their decision making logic. Each controller implements
the following interface:

`predict(state) -> Tuple[action]` or `predict(state) -> action`


"""

from typing import Union
from typing import Union, Tuple
from collections import deque

from sklearn.base import BaseEstimator
@@ -171,17 +178,19 @@ class FeedbackController(BaseEstimator):

class SimpleFeedbackController(BaseEstimator):

    def __init__(self, bounds, stepsize:float=1, window: int=1):
        self.bounds = np.asarray(bounds)
    def __init__(self, bounds, stepsize:float=1, window: int=1, seed=None):
        self.bounds = np.asarray(bounds) # 1D array of (min, max) for setpoint
        self.stepsize = stepsize
        self.window = window
        self.seed = seed
        self.random = np.random.RandomState(seed) # pylint: disable=no-member
        self._feedbacks = deque(maxlen=100)
        self._states = deque(maxlen=100)
        self._actions = deque(maxlen=100)
        self._errors = deque(maxlen=100)


    def predict(self, X: pd.DataFrame):
    def predict(self, X: pd.DataFrame) -> Tuple[np.ndarray]:
        feedback = self.feedback(X)
        self._feedbacks.append(feedback)

@@ -195,21 +204,27 @@ class SimpleFeedbackController(BaseEstimator):
            a_2, a_1 = self._actions[-1], self._actions[-2]
            f_2, f_1 = self._feedbacks[-1], self._feedbacks[-2]
            # What was the direction of change in action from the last 2 steps?
            dir_a = np.random.choice([-1, 1]) if np.sign(a_2 - a_1) == 0 else np.sign(a_2 - a_1)
            dir_a = np.sign(a_2 - a_1)
            # What was the direction of change in feedback from the last 2 steps?
            dir_f = np.sign(f_2 - f_1)
            # Feedback dir, action dir, step action
            #       -           -           +
            #       -           +           -
            #       +           -           -
            #       +           +           +
            #       0           0          rnd
            #       0           -          rnd
            #       0           +          rnd
            #       -           0          rnd
            #       -           -           +       dir_f * dir_a
            #       -           +           -       dir_f * dir_a
            #       +           0           0       dir_f * dir_a
            #       +           -           -       dir_f * dir_a
            #       +           +           +       dir_f * dir_a
            if (dir_f == 0 or (dir_f < 0 and dir_a == 0)):
                step_action = self.stepsize * np.random.choice([-1, 1], size=1)
            else:
                step_action = self.stepsize * dir_a * dir_f
            action = a_2 + step_action
            action = self.clip_action(action, X)
            
            self._actions.append(action)
        # print('err: {:8.2f}, d_err: {:8.2f}, T: {:5.2f}, deltaT: {:5.2f}'\
        #     .format(error, delta_error, action[0], step_action[0]))
        return action,


+1 −1
Changes for src/controller.py: 1 added line, 1 removed line.
Original line number Diff line number Diff line
@@ -197,7 +197,7 @@ def get_controller(**settings) -> BaseEstimator:
                return - X['PowChi'] - X['PowFanA'] - X['PowFanB'] - X['PowConP']
        
        def starting_action(self, X):
            return np.asarray([X['TempWetBulb'] + np.random.uniform(low=4, high=6)])
            return np.asarray([X['TempWetBulb'] + self.random.uniform(low=4, high=6)])

        def clip_action(self, u, X):
            u = super().clip_action(u, X)
Loading