C4_Sweep_Test.py

# -*- coding: utf-8 -*-
"""
Created on Mon Nov 10 12:06:45 2025

@author: danap
"""

import os
import shutil
import datetime as dt
import pandas as pd
from ochre import Dwelling
from ochre.utils.schedule import ALL_SCHEDULE_NAMES
import concurrent.futures
import random
import time
import datetime
import numpy as np


print(datetime.datetime.fromtimestamp(time.time(), datetime.timezone.utc).astimezone().strftime('%Y-%m-%d %H:%M:%S %Z')) 
start_time = time.time()

#########################################
# USER SETTINGS
#########################################

filename = '180113_1_3_BaselineLoad' # date that's thrown away, num of simulation days, data res, ramp or no ramp control
# 04 / 07 


# Paths
DEFAULT_INPUT = r"C:\Users\danap\anaconda3\Lib\site-packages\ochre\defaults\Input Files"
DEFAULT_WEATHER = r"C:\Users\danap\anaconda3\Lib\site-packages\ochre\defaults\Weather\USA_OR_Portland.Intl.AP.726980_TMY3.epw"
WORKING_DIR = r"C:\Users\danap\OCHRE_Working"
INPUT_DIR = os.path.join(WORKING_DIR, "Input Files")
WEATHER_DIR = os.path.join(WORKING_DIR, "Weather")
WEATHER_FILE = os.path.join(WEATHER_DIR, "USA_OR_Portland.Intl.AP.726980_TMY3.epw")

# Simulation parameters
Start = dt.datetime(2018, 1, 13, 0, 0)
Duration = 2  # days
t_res = 3  # minutes
jitter_min = 5

# HPWH control parameters (°F)
Tcontrol_SHEDF = 122
step = 5
ER_list = []
Tcontrol_dbF = np.arange(7, 7 + step, step)  # Deadband sweep list (°F)
Tcontrol_LOADF = 130
Tcontrol_LOADdeadbandF = 2
TbaselineF = 130
TdeadbandF = 7
Tinit = 128

# Base schedule template
my_schedule = {
    'M_LU_time': '03:00',
    'M_LU_duration': 3,
    'M_S_time': '06:00',
    'M_S_duration': 4,
    'E_ALU_time': '14:00',
    'E_ALU_duration': 3,
    'E_S_time': '17:00',
    'E_S_duration': 3
}

# Randomization bins
M_LU_weights = [14, 28, 34, 41, 46, 46, 41, 33, 30, 31, 35, 30]
M_LU_bins = pd.date_range("03:00", periods=len(M_LU_weights), freq="15min").strftime("%H:%M").tolist()
E_ALU_weights = [17, 21, 27, 37, 40, 46, 40, 42, 36, 32, 33, 38]
E_ALU_bins = pd.date_range("14:00", periods=len(E_ALU_weights), freq="15min").strftime("%H:%M").tolist()

#########################################
# TEMPERATURE CONVERSIONS F to C
#########################################

def f_to_c(temp_f): 
    return (temp_f - 32) * 5/9

def f_to_c_DB(temp_f):
    return 5/9 * temp_f

Tcontrol_SHEDC = f_to_c(Tcontrol_SHEDF)
Tcontrol_deadbandC_list = Tcontrol_dbF * 5/9  
Tcontrol_LOADC = f_to_c(Tcontrol_LOADF)
Tcontrol_LOADdeadbandC = f_to_c_DB(Tcontrol_LOADdeadbandF)
TbaselineC = f_to_c(TbaselineF)
TdeadbandC = f_to_c_DB(TdeadbandF)
TinitC = f_to_c(Tinit)

#########################################
# HPWH CONTROL FUNCTION
#########################################

def determine_hpwh_control(sim_time, current_temp_c, sched_cfg, deadband_C, **kwargs):  
    ctrl_signal = {
        'Water Heating': {
            'Setpoint': TbaselineC,
            'Deadband': TdeadbandC,
            'Load Fraction': 1,
        }
    }

    base_date = sim_time.date()
    def get_time_range(key_prefix):
        start = pd.to_datetime(f"{base_date} {sched_cfg[f'{key_prefix}_time']}")
        end = start + pd.Timedelta(hours=sched_cfg[f'{key_prefix}_duration'])
        return start, end

    ranges = {
        'M_LU': get_time_range('M_LU'),
        'M_S': get_time_range('M_S'),
        'E_ALU': get_time_range('E_ALU'),
        'E_S': get_time_range('E_S'),
    }

    if ranges['M_LU'][0] <= sim_time < ranges['M_LU'][1] or ranges['E_ALU'][0] <= sim_time < ranges['E_ALU'][1]:
        ctrl_signal['Water Heating'].update({
            'Setpoint': Tcontrol_LOADC,
            'Deadband': Tcontrol_LOADdeadbandC
        })
    elif ranges['M_S'][0] <= sim_time < ranges['M_S'][1] or ranges['E_S'][0] <= sim_time < ranges['E_S'][1]:
        ctrl_signal['Water Heating'].update({
            'Setpoint': Tcontrol_SHEDC,
            'Deadband': deadband_C  # <<< use the current sweep value
        })

    return ctrl_signal

#########################################
# SCHEDULE FILTERING
#########################################

def filter_schedules(home_path):
    orig_sched_file = os.path.join(home_path, 'schedules.csv')
    filtered_sched_file = os.path.join(home_path, 'filtered_schedules.csv')

    df_sched = pd.read_csv(orig_sched_file)
    valid_schedule_names = set(ALL_SCHEDULE_NAMES.keys())

    hpwh_cols = ['M_LU_time','M_LU_duration','M_S_time','M_S_duration',
                 'E_ALU_time','E_ALU_duration','E_S_time','E_S_duration']
    filtered_columns = [col for col in df_sched.columns if col in valid_schedule_names or col in hpwh_cols]

    dropped_columns = [col for col in df_sched.columns if col not in filtered_columns]
    if dropped_columns:
        print(f"Dropped invalid schedules for {home_path}: {dropped_columns}")

    df_sched_filtered = df_sched[filtered_columns]
    df_sched_filtered.to_csv(filtered_sched_file, index=False)
    return filtered_sched_file

#########################################
# SIMULATION FUNCTION
#########################################

def simulate_home(home_path, weather_file_path, schedule_cfg, deadband_C):  # <<< CHANGED
    filtered_sched_file = filter_schedules(home_path)
    hpxml_file = os.path.join(home_path, 'in.XML')
    results_dir = os.path.join(home_path, "Results")
    os.makedirs(results_dir, exist_ok=True)

    dwelling_args_local = {
        "start_time": Start,
        "time_res": dt.timedelta(minutes=t_res),
        "duration": dt.timedelta(days=Duration),
        "hpxml_file": hpxml_file,
        "hpxml_schedule_file": filtered_sched_file,
        "weather_file": weather_file_path,
        "verbosity": 7,
        "Equipment": {
            "Water Heating": {
                "Initial Temperature (C)": TinitC, 
                "hp_only_mode": True,
                "Max Tank Temperature": 70,
                "Upper Node": 3,
                "Lower Node": 10,
                "Upper Node Weight": 0.75,
            },
        }
    }

    sim_dwelling = Dwelling(name="HPWH Controlled", **dwelling_args_local)
    hpwh_unit = sim_dwelling.get_equipment_by_end_use('Water Heating')

    for sim_time in sim_dwelling.sim_times:
        
        # --- Day 1: force baseline ---
        if sim_time < Start + pd.Timedelta(days=1):
            control_cmd = {
                'Water Heating': {
                    'Setpoint': TbaselineC,
                    'Deadband': TdeadbandC,
                    'Load Fraction': 1,
                }
            }
            sim_dwelling.update(control_signal=control_cmd)
            continue
        
# =============================================================================
#         this is added delete until next block
# =============================================================================
        else:
            control_cmd = {
                'Water Heating': {
                    'Setpoint': TbaselineC,
                    'Deadband': TdeadbandC,
                    'Load Fraction': 1,
                }
            }
            sim_dwelling.update(control_signal=control_cmd)
            continue
        
        

# =============================================================================
#         
#         current_setpt = hpwh_unit.schedule.loc[sim_time, 'Water Heating Setpoint (C)']
#         control_cmd = determine_hpwh_control(sim_time=sim_time,
#                                              current_temp_c=current_setpt,
#                                              sched_cfg=schedule_cfg,
#                                              deadband_C=deadband_C)
#         
#         
#         
# =============================================================================
        
        
        sim_dwelling.update(control_signal=control_cmd)
    df_ctrl, _, _ = sim_dwelling.finalize()

    df_ctrl = remove_first_day(df_ctrl, Start)

    CTRL_COLS = ["Time", "Total Electric Power (kW)",
                 "Total Electric Energy (kWh)",
                 "Water Heating Electric Power (kW)",
                 "Water Heating COP (-)",
                 "Water Heating Deadband Upper Limit (C)",
                 "Water Heating Deadband Lower Limit (C)",
                 "Water Heating Heat Pump COP (-)",
                 "Water Heating Control Temperature (C)",
                 "Hot Water Outlet Temperature (C)",
                 "Temperature - Indoor (C)"]

    df_ctrl = df_ctrl[[c for c in CTRL_COLS if c in df_ctrl.columns]]

    suffix = f"_DB{round(deadband_C * 9/5)}F"  # <<< ADDED: file suffix
    out_path = os.path.join(results_dir, f"hpwh_controlled{suffix}.parquet")
    df_ctrl.to_parquet(out_path, index=False)

    # <<< CHANGED >>> Do NOT delete other sweep results. Keep all _DB*.parquet files.
    # If you want to remove non-parquet junk, do that safely here (example below).
    for item in os.listdir(results_dir):
        path = os.path.join(results_dir, item)
        if os.path.isfile(path) and not item.endswith(".parquet"):
            try:
                os.remove(path)
            except Exception as e:
                print(f"Could not delete {path}: {e}")

    return df_ctrl

#########################################
# FIND ALL HOMES
#########################################

def find_all_homes(base_dir):
    homes = []
    for item in os.listdir(base_dir):
        home_path = os.path.join(base_dir, item)
        if os.path.isdir(home_path):
            if os.path.isfile(os.path.join(home_path, 'in.XML')) and \
               os.path.isfile(os.path.join(home_path, 'schedules.csv')):
                homes.append(home_path)
    return homes

#########################################
# DELETE FIRST DAY ONLY
#########################################

def remove_first_day(df, start_date):
    if 'Time' not in df.columns:
        df = df.reset_index()
        if 'index' in df.columns:
            df.rename(columns={'index': 'Time'}, inplace=True)
    df['Time'] = pd.to_datetime(df['Time'], errors='coerce')
    first_day_end = start_date + pd.Timedelta(days=1)
    return df[df['Time'] >= first_day_end].copy()

#########################################
# CLEAN UP FILES
#########################################

def cleanup_results_dir(results_dir, keep_files=None):
    if keep_files is None:
        keep_files = []
    for item in os.listdir(results_dir):
        path = os.path.join(results_dir, item)
        if os.path.isfile(path) and item not in keep_files:
            try:
                os.remove(path)
            except Exception as e:
                print(f"Could not delete {path}: {e}")
        elif os.path.isdir(path):
            try:
                shutil.rmtree(path)
            except Exception as e:
                print(f"Could not delete folder {path}: {e}")

#########################################
# MAIN EXECUTION
#########################################

if __name__ == "__main__":
    os.makedirs(INPUT_DIR, exist_ok=True)
    os.makedirs(WEATHER_DIR, exist_ok=True)

    for item in os.listdir(DEFAULT_INPUT):
        src = os.path.join(DEFAULT_INPUT, item)
        dst = os.path.join(INPUT_DIR, item)
        if os.path.isdir(src) and not os.path.exists(dst):
            shutil.copytree(src, dst)

    if not os.path.exists(WEATHER_FILE):
        shutil.copy(DEFAULT_WEATHER, WEATHER_FILE)

    homes = find_all_homes(INPUT_DIR)
    print(f"Found {len(homes)} homes")

    # Weighted pools setup (same as before)
    M_LU_weighted_pool = [bin_time for bin_time, weight in zip(M_LU_bins, M_LU_weights) for _ in range(weight)]
    # random.shuffle(M_LU_weighted_pool)

    MS_bins = pd.date_range("10:00", "13:45", freq="15min")
    MS_weights = [20, 23, 24, 23, 22, 22, 25, 26, 26, 29, 29, 29, 29, 27, 28, 27]
    MS_offsets = [(t - pd.Timestamp("10:00")).total_seconds()/3600 for t in MS_bins]
    MS_weighted_pool = [offset for offset, w in zip(MS_offsets, MS_weights) for _ in range(w)]
    # random.shuffle(MS_weighted_pool)

    E_ALU_weighted_pool = [bin_time for bin_time, weight in zip(E_ALU_bins, E_ALU_weights) for _ in range(weight)]
    # random.shuffle(E_ALU_weighted_pool)
    
    ES_bins = pd.date_range("20:00", "23:45", freq="15min")
    ES_weights = [17, 21, 24, 25, 26, 24, 24, 23, 23, 23, 23, 25, 28, 30, 33, 40]
    ES_offsets = [(t - pd.Timestamp("20:00")).total_seconds()/3600 for t in ES_bins]
    ES_weighted_pool = [offset2 for offset2, m in zip(ES_offsets, ES_weights) for _ in range(m)]
    # random.shuffle(ES_weighted_pool)

    # Assign schedules per home (unchanged)
    home_schedules = {}
    fmt = "%H:%M"
    for home in homes:
        sched = my_schedule.copy()
        if M_LU_weighted_pool:
            M_LU_base = M_LU_weighted_pool.pop()
        else:
            M_LU_base = random.choice(M_LU_bins)
        t_base = pd.to_datetime(M_LU_base, format=fmt)
        jitter = pd.Timedelta(minutes=random.uniform(-jitter_min, jitter_min))
        t_jittered = t_base + jitter
        sched['M_LU_time'] = t_jittered.strftime(fmt)

        if M_LU_base == '05:45':
            t_MS_start = pd.to_datetime("06:15", format=fmt)
        else:
            t_MS_start = pd.to_datetime(my_schedule['M_S_time'], format=fmt)
        t_MS_start += pd.Timedelta(minutes=random.uniform(-jitter_min, jitter_min))
        sched['M_S_time'] = t_MS_start.strftime(fmt)

        t_MLU_start = pd.to_datetime(sched['M_LU_time'], format=fmt)
        t_MLU_end = t_MS_start
        if t_MLU_end <= t_MLU_start:
            t_MLU_end += pd.Timedelta(days=1)
        sched['M_LU_duration'] = max(1, (t_MLU_end - t_MLU_start).total_seconds() / 3600)

        if MS_weighted_pool:
            n = MS_weighted_pool.pop()
        else:
            n = random.choice(MS_offsets)
        sched['M_S_duration'] = 4 + n

        if E_ALU_weighted_pool:
            E_ALU_base = E_ALU_weighted_pool.pop()
        else:
            E_ALU_base = random.choice(E_ALU_bins)
        t_E_ALU_start = pd.to_datetime(E_ALU_base, format=fmt)
        jitter = pd.Timedelta(minutes=random.uniform(-jitter_min, jitter_min))
        t_E_ALU_start += jitter
        sched['E_ALU_time'] = t_E_ALU_start.strftime(fmt)

        t_ES_start = pd.to_datetime(my_schedule['E_S_time'], format=fmt)
        t_ES_start += pd.Timedelta(minutes=random.uniform(-jitter_min, jitter_min))
        sched['E_S_time'] = t_ES_start.strftime(fmt)
        
        if t_ES_start <= t_E_ALU_start:
            t_ES_start += pd.Timedelta(days=1)
        sched['E_ALU_duration'] = max(1, (t_ES_start - t_E_ALU_start).total_seconds() / 3600)
        
        if ES_weighted_pool:
            n = ES_weighted_pool.pop()
        else:
            n = random.choice(ES_offsets)
        sched['E_S_duration'] = 3 + n
        home_schedules[home] = sched

    #########################################
    # SWEEP DEADBANDBAND VALUES
    #########################################

    for deadband_C in Tcontrol_deadbandC_list:
        print(f"\n=== Running simulations for deadband {deadband_C:.2f} °C ({round(deadband_C*9/5)} °F) ===")
        with concurrent.futures.ThreadPoolExecutor(max_workers=8) as executor:
            futures = [
                executor.submit(simulate_home, home, WEATHER_FILE, home_schedules[home], deadband_C)
                for home in homes
            ]
            for f in concurrent.futures.as_completed(futures):
                try:
                    f.result()
                except Exception as e:
                    print("Simulation failed:", e)

    #########################################
    # AGGREGATE RESULTS
    #########################################

    def aggregate_results(homes, work_dir):
        for deadband_C in Tcontrol_deadbandC_list:
            suffix = f"_DB{round(deadband_C * 9/5)}F"
            all_ctrl = []
            for home in homes:
                results_dir = os.path.join(home, "Results")
                ctrl_file = os.path.join(results_dir, f"hpwh_controlled{suffix}.parquet")
                if os.path.exists(ctrl_file):
                    df_ctrl = pd.read_parquet(ctrl_file)
                    df_ctrl["Home"] = os.path.basename(home)
                    df_ctrl["Deadband_C"] = deadband_C
                    all_ctrl.append(df_ctrl)
            if all_ctrl:
                df_ctrl_all = pd.concat(all_ctrl, ignore_index=True)
                outp = os.path.join(work_dir, filename + f"{suffix}_Control.parquet")
                df_ctrl_all.to_parquet(outp, index=False)
                print(f"Aggregated results written for {suffix}")

        # # <<< ADDED >>> Combine all per-deadband aggregated files into one master file
        # master_files = [os.path.join(work_dir, f) for f in os.listdir(work_dir)
        #                 if f.endswith("_Control.parquet") and "_DB" in f]
        # if master_files:
        #     all_dfs = [pd.read_parquet(f) for f in master_files]
        #     df_all = pd.concat(all_dfs, ignore_index=True)
        #     master_out = os.path.join(work_dir, filename + "_Control.parquet")
        #     df_all.to_parquet(master_out, index=False)
        #     print(f"Master aggregated file written: {master_out}")

    aggregate_results(homes, WORKING_DIR)
    
    
    # def aggregate_by_setpoint(work_dir):
    #     """
    #     Aggregate all per-deadband HPWH control files for the same setpoint (ShedXXX)
    #     into a single master file per setpoint.
    #     """
    #     # List all files in the directory
    #     all_files = [f for f in os.listdir(work_dir) if f.endswith("_Control.parquet") and "_DB" in f]
    
    #     # Group files by setpoint using regex
    #     setpoint_groups = {}
    #     for f in all_files:
    #         match = re.search(r"(Shed\d+)_DB\d+F_Control\.parquet", f)
    #         if match:
    #             setpoint = match.group(1)
    #             setpoint_groups.setdefault(setpoint, []).append(os.path.join(work_dir, f))
    
    #     # Aggregate each setpoint group
    #     for setpoint, files in setpoint_groups.items():
    #         all_dfs = []
    #         for file_path in files:
    #             df = pd.read_parquet(file_path)
    #             # Optional: keep track of deadband in df
    #             df["SourceFile"] = os.path.basename(file_path)
    #             all_dfs.append(df)
    #         df_master = pd.concat(all_dfs, ignore_index=True)
    #         master_file = os.path.join(work_dir, f"{setpoint}_Control.parquet")
    #         df_master.to_parquet(master_file, index=False)
    #         print(f"Aggregated {len(files)} files for {setpoint} → {master_file}")
            
    # aggregate_by_setpoint(WORKING_DIR)


    def aggregate_current_setpoint(work_dir, prefix):
        """
        Aggregate all per-deadband HPWH control files produced by aggregate_results
        for the same simulation prefix (date, duration, setpoint).
    
        Example:
            prefix = '180110_1_3_Shed110'
            Input files: 180110_1_3_Shed110_DB5F_Control.parquet,
                         180110_1_3_Shed110_DB10F_Control.parquet, ...
            Output file: 180110_1_3_Shed110_Control.parquet
        """
        # Find all _DB*.parquet files for this prefix
        all_files = [
            f for f in os.listdir(work_dir)
            if f.endswith("_Control.parquet") and f.startswith(prefix + "_DB")
        ]
    
        if not all_files:
            print(f"No deadband files found for prefix {prefix}")
            return
    
        all_dfs = []
        for f in all_files:
            path = os.path.join(work_dir, f)
            df = pd.read_parquet(path)
            df["SourceFile"] = f  # track which deadband this came from
            all_dfs.append(df)
    
        df_master = pd.concat(all_dfs, ignore_index=True)
        master_file = os.path.join(work_dir, f"{prefix}_Control.parquet")
        df_master.to_parquet(master_file, index=False)
        print(f"Aggregated {len(all_files)} deadband files for {prefix} → {master_file}")


    aggregate_current_setpoint(WORKING_DIR, filename)


end_time = time.time()
execution_time = end_time - start_time
execution_min = execution_time/60
print(f"Execution time: {execution_min} minutes")