<jupyter_start><jupyter_text>MrCool 3 ton unit ![Screen%20Shot%202022-03-16%20at%208.22.03%20PM.png](attachment:Screen%20Shot%202022-03-16%20at%208.22.03%20PM.png)<jupyter_code>mrcool_specs_3_ton = [     
    [-30, 18400, 4150],
    [-26, 21100, 4430],
    [-21, 25200, 4630],
    [-15, 30000, 4800],
    [-8, 36000, 5450],
    [0, 36000, 4800],
    [8, 38000, 3600],
    [15, 38000, 3150],
    [20, 38000, 3120]
]

mrcool_specs_5_ton = [     
    [-30, 27000, 6350],
    [-26, 32000, 6700],
    [-21, 37000, 7000],
    [-15, 44000, 7160],
    [-8, 49000, 7100],
    [0, 54000, 6600],
    [8, 56000, 5300],
    [15, 60000, 5400],
    [20, 60000, 5100]
]

# numpy.interp(x, xp, fp, left=None, right=None, period=None)

btuph_to_w = 0.000293071 * 1000.

mrcool_specs = mrcool_specs_3_ton

def mrcool_heat_output(t_out):
    h_out_btuph = np.interp(t_out, [x[0] for x in mrcool_specs], [x[1] for x in mrcool_specs])
    h_out_w = btuph_to_w * h_out_btuph
    return round(h_out_w, 1)

def mrcool_power_consumption(t_out):
    power_consumption_w = np.interp(t_out, [x[0] for x in mrcool_specs], [x[2] for x in mrcool_specs])
    return power_consumption_w

def electricity_cost(hour_of_day, day_of_week):
    # https://www.hydroone.com/rates-and-billing/rates-and-charges/electricity-pricing-and-costs
    
    off_peak = 8.2
    mid_peak = 11.3
    on_peak = 17.0
    
    if day_of_week in (5, 6):
        # weekend
        return off_peak
    
    if hour_of_day < 7 or hour_of_day > 19:
        return off_peak
    
    if hour_of_day > 11 and hour_of_day < 17:
        return mid_peak
    
    return on_peak
kwargs = {"header": 4, "index_col": False}
ecobee = pd.concat([
    pd.read_csv("ecobee_data/2021-09.csv", **kwargs),
    pd.read_csv("ecobee_data/2021-10.csv", **kwargs),
    pd.read_csv("ecobee_data/2021-11.csv", **kwargs),
    pd.read_csv("ecobee_data/2021-12.csv", **kwargs),
    pd.read_csv("ecobee_data/2022-01.csv", **kwargs),
    pd.read_csv("ecobee_data/2022-02.csv", **kwargs),
    pd.read_csv("ecobee_data/2022-03-partial.csv", **kwargs)
])

ecobee.columns
ecobee["datetime"] = pd.to_datetime(ecobee["Date"] + "T" + ecobee["Time"])
cols_of_interest = [
    "datetime",
    "System Setting",
    "System Mode",
    "Heat Set Temp (C)",
    'Outdoor Temp (C)',
    "Thermostat Temperature (C)",
    'Heat Stage 1 (sec)'
]
col_dict = {
    "System Setting": "system_setting", 
    "System Mode": "system_mode", 
    "Calendar Event": "calendar_event",
    "Outdoor Temp (C)": "outdoor_temp_c",
    "Thermostat Temperature (C)": "indoor_temp_c",
    "Heat Stage 1 (sec)": "heat_1_sec",
    'Heat Set Temp (C)': "heat_setpoint_c"
}


ecobee = ecobee[cols_of_interest].copy()
ecobee.set_index("datetime", inplace=True)
ecobee.rename(columns=col_dict, inplace=True)
ecobee.head()
ecobee["dt_c"] = ecobee["indoor_temp_c"] - ecobee["outdoor_temp_c"]
heating_data = ecobee.query("system_setting == 'heat'").copy()
agg_fns = {"outdoor_temp_c": np.nanmean, "dt_c": np.nanmean, "heat_1_sec": np.nansum, "system_setting": len}
df_h = heating_data.groupby(pd.Grouper(freq="H")).aggregate(agg_fns)
df_h.rename(columns={"system_setting": "data_count"}, inplace=True)
df_h = df_h.query("data_count == 12").copy() # throw out incomplete data
df_h["runtime_fraction"] = df_h["heat_1_sec"]/(60.*60.)
furnace_output_btuph = 62000
df_h["heat_load_w"] = df_h["runtime_fraction"] * furnace_output_btuph * btuph_to_w
x_data = df_h["dt_c"]
y_data = df_h["heat_load_w"]
plt.plot(x_data, y_data, '.')
increment = 0.5
_temp_range = np.arange(0, 50, increment)
temp_range = []

mean_load_w_array = []
for t in _temp_range:
    t1 = t - increment/2.
    t2 = t + increment/2.
    mean_load_w = df_h.query("dt_c >= @t1 and dt_c < @t2")["heat_load_w"].mean()  
    if mean_load_w == mean_load_w:
        # filter out nan, if no t_out exists
        temp_range.append(t1)
        mean_load_w_array.append(mean_load_w)

plt.plot(temp_range, mean_load_w_array, '.')
stats.linregress(temp_range, np.array(mean_load_w_array))

stats.linregress(temp_range, np.array(mean_load_w_array))

slope, intercept, r_value, p_value, std_err = stats.linregress(temp_range, mean_load_w_array)


def get_load_w(dt_c):
    y = slope * dt_c + intercept
    return y

calc_y_data = [get_load_w(x) for x in temp_range]
plt.plot(
    temp_range, mean_load_w_array, '.',
    temp_range, calc_y_data, '-'
)<jupyter_output><empty_output><jupyter_text>Estimate the cost<jupyter_code>agg_fns = {"outdoor_temp_c": np.nanmean, "dt_c": np.nanmean}
df_h_sim = heating_data.groupby(pd.Grouper(freq="H")).aggregate(agg_fns)

df_h_sim = df_h_sim.query("outdoor_temp_c == outdoor_temp_c and dt_c == dt_c").copy()
df_h_sim["modelled_load_w"] = df_h_sim["dt_c"].map(get_load_w)
df_h_sim["heat_pump_max_out_w"] = df_h_sim["outdoor_temp_c"].map(mrcool_heat_output)
df_h_sim["cop"] = df_h_sim["heat_pump_max_out_w"] / df_h_sim["outdoor_temp_c"].map(mrcool_power_consumption)
def out_fn(x):
    if x["heat_pump_max_out_w"] >= x["modelled_load_w"]:
        return x["modelled_load_w"]
    return x["heat_pump_max_out_w"]

df_h_sim["heat_pump_out_w"] = df_h_sim.apply(out_fn, axis=1)
df_h_sim["aux_out_w"] = df_h_sim["modelled_load_w"] - df_h_sim["heat_pump_out_w"]
df_h_sim["heat_pump_power_kwh"] = (1./1000.) * df_h_sim["heat_pump_out_w"] /  df_h_sim["cop"]
df_h_sim["datetime2"] = df_h_sim.index
df_h_sim["electricity_rate"] = df_h_sim["datetime2"].map(lambda x: electricity_cost(x.hour, x.day_of_week))
df_h_sim["heat_pump_cost"] = df_h_sim["heat_pump_power_kwh"] * df_h_sim["electricity_rate"] * (1./100.)
df_h_sim["aux_power_kwh"] = df_h_sim["aux_out_w"] / 1000.
df_h_sim["aux_cost"] = df_h_sim["aux_power_kwh"] * df_h_sim["electricity_rate"] * (1./100.)
aux_cost = df_h_sim["aux_cost"].sum()
aux_cost
heat_pump_cost = df_h_sim["heat_pump_cost"].sum()
heat_pump_cost
print(f"aux cost: \n{round(100 * aux_cost / (aux_cost + heat_pump_cost), 1)}% of total\n")
print(f"heat pump cost: \n{round(100 * heat_pump_cost / (aux_cost + heat_pump_cost), 1)}% of total")
max(df_h_sim["aux_out_w"])
qs = "datetime2 < '2022-02-01'"
df_h_sim.query(qs)["aux_cost"].sum() + df_h_sim.query(qs)["heat_pump_cost"].sum()
df_h_sim.head()<jupyter_output><empty_output><jupyter_text>Load Enbridge data<jupyter_code>enbridge = pd.read_csv("./enbridge_data.csv")
enbridge.columns
enbridge_col_dict = {
    " Billing Period From": "start_date",
    " Billing Period To": "end_date",
    " Consumption": "gas_m3",
    " Gas Charge": "gas_cost",
    " Invoice Amount": "invoice_amount"
}

enbridge = enbridge[enbridge_col_dict.keys()].rename(columns=enbridge_col_dict).copy()

enbridge["start_date"] = pd.to_datetime(enbridge["start_date"])
enbridge["end_date"] = pd.to_datetime(enbridge["end_date"])
def convert_to_num(x):
    return float(x.strip().replace("$", ""))

enbridge["gas_cost"] = enbridge["gas_cost"].map(convert_to_num)
# enbridge["invoice_amount"] = enbridge["invoice_amount"].map(convert_to_num)
enbridge["invoice_amount"].sum()
enbridge["invoice_amount"].sum()
enbridge<jupyter_output><empty_output><jupyter_text>See what my model predicts for gas usage<jupyter_code>def get_nat_gas_consumption(load_Wh):
    # from Lennox spec sheet
    heat_output_kw = 18.2
    gas_input_kw = 19.3
    efficiency = heat_output_kw / gas_input_kw
    
    # https://www.nrcan.gc.ca/energy/energy-sources-distribution/natural-gas/natural-gas-primer/5641#conversion
    m3_in_GJ = 0.0373
    # 1 kWh = 0.0036 GJ
    # 1m3 nat gas = 0.0373 GJ
    kWh_in_GJ = 0.0036 # google
    kWh_in_m3 = kWh_in_GJ / m3_in_GJ
    
    load_kWh = load_Wh / 1000.
    input_kWh = load_kWh * (1./ efficiency)
    input_m3 = input_kWh * kWh_in_m3
    return input_m3
df_h_sim["ref_gas_cons_m3"] = df_h_sim["modelled_load_w"].map(get_nat_gas_consumption)
df_h_sim.head()
def get_gas_consumed(start_date, end_date):
    
    return df_h_sim["ref_gas_cons_m3"]["2021-09-29":"2021-10-01"].sum()
def get_modelled_gas_cons(x):
    start_date = x["start_date"]
    end_date = x["end_date"]
    return df_h_sim["ref_gas_cons_m3"][start_date:end_date].sum()    

enbridge["modelled_gas_m3"] = enbridge.apply(get_modelled_gas_cons, axis=1)
enbridge
enbridge["modelled_pct_error"] = (enbridge["gas_m3"] - enbridge["modelled_gas_m3"]) / enbridge["gas_m3"]
season_start_date = '2021-10-23'
season_end_date = '2022-02-22'

heating_season = enbridge.query("start_date >= @season_start_date").copy()
heating_season["gas_cost"].sum()
sim_total = df_h_sim["heat_pump_cost"][season_start_date:season_end_date].sum() + df_h_sim["aux_cost"][season_start_date:season_end_date].sum()
sim_total
correction = (heating_season["gas_m3"].sum() - heating_season["modelled_gas_m3"].sum()) / heating_season["modelled_gas_m3"].sum()

correction
corr_sim_total = sim_total * (1 + correction)
corr_sim_total
(corr_sim_total - heating_season["gas_cost"].sum() )/ heating_season["gas_cost"].sum()
mrcool_heat_output(-15)/mrcool_power_consumption(-15)<jupyter_output><empty_output><jupyter_text>Greenhouse gas emmissions<jupyter_code>total_hp_kwh = df_h_sim["heat_pump_power_kwh"][season_start_date:season_end_date].sum() + df_h_sim["aux_power_kwh"][season_start_date:season_end_date].sum()
modelled_gas_m3 = heating_season["modelled_gas_m3"].sum()
furnace_tco2 = 0.002 * modelled_gas_m3
print(furnace_tco2)
hp_tco2 = 0.001*0.24*total_hp_kwh
print(hp_tco2)
(furnace_tco2 - hp_tco2) / furnace_tco2
df_h_sim.head()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15661881977671452, 'fraction_numerical': 0.0531632110579479, 'mean_word_length': 5.304289544235925, 'pattern_counts': {'":': 21, '<': 25, '<?xml version=': 0, '>': 26, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '447031', 'license': 'MIT', 'n_tokens_mistral': 4381, 'n_tokens_neox': 4002, 'n_words': 578, 'path': 'heat_pump.ipynb', 'repo_name': 'smcmurtry/heat_pump', 'size': '9405'}