<jupyter_start><jupyter_text>Analysis of 2015 CCHRC Air Source Heat Pump Study Data<jupyter_code>from datetime import datetime, timedelta
import pandas as pd
import numpy as np
# import matplotlib pyplot commands
from matplotlib.pyplot import *
# Show Plots in the Notebook
%matplotlib inline
 
# 'style' the plot like fivethirtyeight.com website
style.use('bmh')
rcParams['figure.figsize']= (9, 6)   # set Chart Size
rcParams['font.size'] = 14            # set Font size in Chart<jupyter_output><empty_output><jupyter_text>Dillingham Mitsubishi Mini-Split Analysis<jupyter_code>fn = 'data/Dillingham_DLG_DHP_TenSecond.dat'
# Take a look at the header rows
f = open(fn)
for i in range(5):
    print(f'Line {i}: {f.readline().strip()}')
dfd = pd.read_csv(fn, skiprows=[0, 2, 3], parse_dates=['TIMESTAMP'], na_values=['NAN'])
dfd.head()
dfd.info()
# Look at some statistics for periods when the unit is producing heat
# Shows the differences between the 3 Supply and 3 Return temperature sensors.
dfd[dfd.DLG_DeliveryT2_F > 90].describe()
# Calculate average delivery and return temperatures, and clean up DataFrame
# renaming columns to shorter names.
dfd['t_delivery'] = (dfd.DLG_DeliveryT1_F + dfd.DLG_DeliveryT2_F + dfd.DLG_DeliveryT3_F) / 3.0
dfd['t_return'] = (dfd.DLG_ReturnT1_F + dfd.DLG_ReturnT2_F + dfd.DLG_ReturnT3_F) / 3.0
dfd.drop(columns=['RECORD', 'DLG_DeliveryT1_F', 'DLG_DeliveryT2_F', 'DLG_DeliveryT3_F', 'DLG_ReturnT1_F', 'DLG_ReturnT2_F', 'DLG_ReturnT3_F'], inplace=True)
dfd.rename(columns={
    'DLG_CoilT_F': 't_coil',
    'DLG_OutsideT_F': 't_out',
    'DLG_Electricity_W_Avg': 'power',
    'DLG_Airflow_fpm_Avg': 'airspeed',
    'TIMESTAMP': 'ts',
}, inplace=True)
dfd.set_index('ts', inplace=True)
dfd.describe()
# Check calibration of Delivery and Return Temp Sensors.
# They should be pretty close to each other when heat pump is off.
df_off = dfd.query('power < 100')
delta_t = df_off.t_delivery - df_off.t_return
# also knock out large delta-T's
delta_t = delta_t.where(np.abs(delta_t) < 3)
print(f'Average Delivery delta-T when Off: {delta_t.mean():.2f}')
delta_t.hist(bins=30)
# How important is accurate delta-T?  What is the average delta-T
# when the unit is On?
df_on = dfd.query('power >= 100')
delta_t = df_on.t_delivery - df_on.t_return
print(delta_t.mean())

# About 38 deg delta-T.  So inaccuracy of 0.3 out of 38 F is about
# 1%.  No big deal if you are looking at the COP just during the On
# periods.
# Compute Air Flow in CFM from calibration data
# Calibration data from:  data/scripts/DLG - calibrating proxy 2015-01-11 - final.xlsx
airflow = 1.3083 * dfd.airspeed + 5.1391         # Air flow in cfm from proxy calibration
nonzero_speed = (dfd.airspeed > 40)              # Speed less than 40 fpm counts as zero speed
dfd['airflow'] = airflow * nonzero_speed         # Zero speed means zero flow
dfd.airflow.hist(bins=30)
# Look like there may be some outliers
dfd.query('airflow > 250').airflow.hist(bins=30)
# Substitute NaN for values over 330
dfd['airflow'] = dfd.airflow.where(dfd.airflow < 330)
dfd.query('airflow > 0').airflow.hist(bins=30)
xlabel('Air Flow, cfm')
ylabel('# of Samples')
axvline(x=286, linestyle='--', color='red')
# Very few points reach the high speed calibration value.  Unit ran on
# High a fair amount, so should see value in that vicinity.
# Also, Mitsubishi Specs for Air Flow are:
#   Powerful: 381 cfm
#   High    : 367 cfm
#   Medium  : 240 cfm
#   Low     : 166 cfm
# Quite a bit higher than the peak observed in the following graph.
# Although other researchers have observed flow rates less than 
# specs as well.
# Calculate Heat Output in Watts

# Constants:
c_p = 0.240 # Specific heat of air in Btu/lb-F
ro = 0.0749   # Density of air (in lb/ft^3) at standard conditions
VHC = c_p*ro  # Volumetric heat capacity of air in Btu/ft^3-F

dfd['heat_out'] = dfd.airflow * (dfd.t_delivery - dfd.t_return) * VHC * 60 / 3.412

# Eliminate the coil temp column
dfd2 = dfd.drop(columns=['t_coil'])

# drop the rows that have one or more NaNs
dfd2.dropna(inplace=True)

dfd2.describe()
# Looking at the overall mean power and heat_out, COP is 1.56 for the
# entire period, counting periods when the heat pump is Off.
# That matches the value on page 48 of the report.
# Date to start for the graph in the next cell.
dt = datetime(2014, 10, 5)
# Step through the days looking at the heat pump power consumption
# across the day.
dfd[dt.strftime('%Y-%m-%d')].power.plot()
dt += timedelta(days=1)
dfd2.query('power > 50').power.hist(bins=30)
# Look at Power Consumption when the heat pump is essentially off
dfd2.query('power <= 50').power.hist(bins=60)
dfd2.query('airspeed > 40').head()
dfd2.describe()
# Function to average over some 'interval_string' time period,
# and then plot Power against outdoor temperature.
def plot_avg_power(df, interval_string):
    df_avg = df.resample(interval_string).mean()
    df_avg.plot(x='t_out', y='power', kind='scatter')
# Experiment with different averaging intervals to see where the data settles down.
plot_avg_power(dfd2, '2D')
# Do 1 day averages and compute COP for those intervals
# Look at filtering out the low power readings so that the COP
# only considers time periods when the heat pump is running. Looking
# at On periods is more consistent with operation in a normal load home,
# as opposed to the very low load home Tom has.
# NEED TO NOTE this when I present this graph.
power_threshold = 25   # Watts required to consider the unit running
dfd2_1day = dfd2.query('power >= @power_threshold').resample('1D').mean()
dfd2_1day['cop'] = dfd2_1day.heat_out / dfd2_1day.power
dfd2_1day.plot(x='t_out', y='cop', kind='scatter')
xlabel('Outdoor Temperature, ℉')
ylim(0,3)
# Calculate seasonal average COP, only looking at periods when heat
# pump is on.
dfd2_mean = dfd2.query('power >= @power_threshold').mean()
dfd2_mean.heat_out / dfd2_mean.power<jupyter_output><empty_output><jupyter_text>Create Standard COP vs. Binned Temperature Plot<jupyter_code>dfb = dfd2_1day.copy()

# calculate the bin for each row.  4 degree bins
dfb['temp_bin'] = dfb['t_out'] // 4.0

# average by bins
dfa = dfb.groupby('temp_bin').mean()
dfa = dfa[['t_out', 'cop']]
dfa.columns = ['temperature', 'cop']
plot(dfa.temperature, dfa.cop, 'b-', marker='o', markersize=5)
plot(dfd2_1day.t_out, dfd2_1day.cop, linestyle='None', color='#555555', marker='+', markersize=7)
xlabel('Outdoor Temperature, ℉')
ylabel('COP')
ylim(0, 3.0)

dfa['source'] = 'CCHRC - Dillingham'
dfa['field_study'] = True
dfa['hspf'] = 10.0
dfa.to_pickle('temp_cop_cchrc_dillingham.pkl')
dfa<jupyter_output><empty_output><jupyter_text>Wrangell Fujitsu Mini-Split Analysis<jupyter_code>fn = 'data/Wrangell_WRG_DHP_TenSecond.dat'
dfd = pd.read_csv(fn, skiprows=[0, 2, 3], parse_dates=['TIMESTAMP'], na_values=['NAN'])
dfd.head()
dfd.info()
# Calculate average delivery and return temperatures, and clean up DataFrame
# renaming columns to shorter names.
dfd['t_delivery'] = (dfd.WRG_DeliveryT1_F + dfd.WRG_DeliveryT2_F + dfd.WRG_DeliveryT3_F) / 3.0
dfd['t_return'] = (dfd.WRG_ReturnT1_F + dfd.WRG_ReturnT2_F + dfd.WRG_ReturnT3_F) / 3.0
dfd.drop(columns=['RECORD', 'WRG_DeliveryT1_F', 'WRG_DeliveryT2_F', 'WRG_DeliveryT3_F', 'WRG_ReturnT1_F', 'WRG_ReturnT2_F', 'WRG_ReturnT3_F'], inplace=True)
dfd.rename(columns={
    'WRG_CoilT_F': 't_coil',
    'WRG_OutsideT_F': 't_out',
    'WRG_Electricity_W_Avg': 'power',
    'WRG_Airflow_fpm_Avg': 'airspeed',
    'TIMESTAMP': 'ts',
}, inplace=True)
dfd.set_index('ts', inplace=True)
# Report states that the valid data period is 12/3/2014 - 5/7/2015
dfd = dfd['2014-12-03' : '2015-05-07'].copy()
dfd.head()
dfd.describe()
# Compute Air Flow in CFM from calibration data
# Calibration data from:  data/scripts/ASHP_WRG.m
airflow = 1.1922 * dfd.airspeed - 26.052         # Air flow in cfm from proxy calibration
nonzero_speed = (dfd.airspeed > 40)              # Speed less than 40 fpm counts as zero speed
dfd['airflow'] = airflow * nonzero_speed         # Zero speed means zero flow
dfd.airflow.hist(bins=30)
# Calculate Heat Output in Watts

# Constants:
c_p = 0.240 # Specific heat of air in Btu/lb-F
ro = 0.0749   # Density of air (in lb/ft^3) at standard conditions
VHC = c_p*ro  # Volumetric heat capacity of air in Btu/ft^3-F

dfd['heat_out'] = dfd.airflow * (dfd.t_delivery - dfd.t_return) * VHC * 60 / 3.412

# Eliminate the coil temp column
dfd2 = dfd.drop(columns=['t_coil'])

# drop the rows that have one or more NaNs
dfd2.dropna(inplace=True)

dfd2.describe()
# Looking at the overall mean power and heat_out, COP is 4.52 for the
# entire period.  That matches the value on page 48 of the report.
dfd2.power.hist(bins=30)
# Experiment with different averaging intervals to see where the data settles down.
plot_avg_power(dfd2, '1D')
ylim(0)
# Do 1 day averages and compute COP for those intervals
dfd2_1day = dfd2.resample('1D').mean()
dfd2_1day['cop'] = dfd2_1day.heat_out / dfd2_1day.power
dfd2_1day.plot(x='t_out', y='cop', kind='scatter')
xlabel('Outdoor Temperature, ℉')
# Calculate seasonal average COP
# This matches the value in the report on page 48.
dfd2_mean = dfd2.mean()
dfd2_mean.heat_out / dfd2_mean.power<jupyter_output><empty_output><jupyter_text>Create Standard COP vs. Binned Temperature Plot<jupyter_code>dfb = dfd2_1day.copy()

# calculate the bin for each row.  4 degree bins
dfb['temp_bin'] = dfb['t_out'] // 4.0

# average by bins
dfbin = dfb.groupby('temp_bin').mean()
dfa = dfbin[['t_out', 'cop']]
dfa.columns = ['temperature', 'cop']
plot(dfa.temperature, dfa.cop, 'b-', marker='o', markersize=5)
plot(dfd2_1day.t_out, dfd2_1day.cop, linestyle='None', color='#555555', marker='+', markersize=7)
xlabel('Outdoor Temperature, ℉')
ylabel('COP')<jupyter_output><empty_output><jupyter_text>The downward slope of COP above 46 F may be due to cycling.  The report shows an example of cycling on page 46.  That sample occurred at 49.4 F outdoors, a temperature in the downward portion of the above curve.  The Kaluza heat pump (without the external thermostat) showed the same characteristic, but a lower outdoor temperature due to the low load of the building.<jupyter_code>dfa['source'] = 'CCHRC - Wrangell'
dfa['field_study'] = True
dfa['hspf'] = 12.0
dfa.to_pickle('temp_cop_cchrc_wrangell.pkl')
dfa<jupyter_output>/home/tabb99/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:1: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  """Entry point for launching an IPython kernel.
/home/tabb99/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:2: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
  
/home/tabb99/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the docum[...]<jupyter_text>Look at the Amount of Heat Load Accessed by the Unit This is done in order to estimate the amount of heat loadaccessed by the unit.  An important issue in the heat pump modelingis how the heat distribution of the space limits the amount of heatthat can be used from the heat pump.<jupyter_code>plot(dfd2_1day.t_out, dfd2_1day.heat_out, linestyle='None', color='#555555', marker='+', markersize=7)
xlabel('Outdoor Temperature, ℉')
ylabel('Heat Output, Watts')
plot(dfbin.t_out, dfbin.heat_out, 'b-', marker='o', markersize=5)
ylim(0)
# Look at the slope of those various segments to calculate the
# the incremental heat loss coefficient of the space served
# by the heat pump.
dfbin['ua'] = - dfbin.heat_out.diff() / dfbin.t_out.diff()
dfbin[['t_out', 'ua']]
# Not very conclusive.  At warm temperature (46 - 53 F), looks like space has
# a pretty large UA:  112 W / F or 382 Btu/hour/F.
# Ignoring a number of the low output points that occur from 20 - 35 F,
# there seems to be a slope going from (27 F, 3000 W) to (52 F, 700 W).
# This is a UA of:
ua = (3000 - 700) / (52 - 27)
f'{ua:.0f} W/F, {ua * 3.412:.0f} Btu/hr/F'
# Instead, lets look at the Heat Output at 25 F, by interpolating the
# binned values.
out_25 = np.interp(25.0, dfbin.t_out, dfbin.heat_out)
print(f'{out_25:.0f} W, {out_25 * 3.412:.0f} Btu/hr')

# Also show the implied UA.  The balance point for this building looks
# like 60 F, so delta-T is 35 F.
print(f'{out_25 * 3.413 / 35.0:.0f} Btu/hr/F')

# We are not doing this for the Dillingham building because that
# building is very low load, and the output at 25 F is surely the
# total building load, not limited by heat distribution issues.<jupyter_output>2164 W, 7384 Btu/hr
211 Btu/hr/F
{'fraction_non_alphanumeric': 0.11936278282284131, 'fraction_numerical': 0.03493920270894259, 'mean_word_length': 4.73477493380406, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '294284', 'license': 'MIT', 'n_tokens_mistral': 5090, 'n_tokens_neox': 4706, 'n_words': 1457, 'path': 'field_studies/cchrc_2015/cchrc_raw_data.ipynb', 'repo_name': 'alanmitchell/heat-pump-study', 'size': '12994'}