# %% 
# Binning of eChiller Power Values
import pandas as pd
import matplotlib.pyplot as plt

# Load the dataset
chiller_data = pd.read_csv('Merge_chiller_15_310823.csv')


P_el_max = 14.5 # kW round 15kW aus Datenblatt 
P_el_min = 0.3 # kW rounded 0 kW aus Datenblatt

Temp_Soll_Netz = 15 # °C Ist Fix Definiert durch anlagenbetreiber definiert

# Define the power range and the number of power states again
power_range = (0, 15)
power_states_count = 5
state_range = 3  # kW

# Binning to compute the intervals for each power state
power_states = {}
for i in range(power_states_count):
    power_states[f"Powerstate {i+1}"] = (i * state_range, (i + 1) * state_range)

print(power_states)
# Define power states and colors
power_states = {
    "Powerstate 1": (0, 3),
    "Powerstate 2": (3, 6),
    "Powerstate 3": (6, 9),
    "Powerstate 4": (9, 12),
    "Powerstate 5": (12, 15)
}
colors = {
    "Powerstate 1": "darkgray",
    "Powerstate 2": "silver",
    "Powerstate 3": "gray",
    "Powerstate 4": "dimgray",
    "Powerstate 5": "lightgray"
}

# Convert the power values to power state labels
power_state_labels = []
for value in chiller_data['WKM_eCH45150Leistung_elek_aktuell']:
    label_assigned = False
    for state, (low, high) in power_states.items():
        if low <= value < high:
            power_state_labels.append(state)
            label_assigned = True
            break
    if not label_assigned:
        power_state_labels.append('Unknown')

# Add the power state labels to the dataframe
chiller_data['Power_State'] = power_state_labels



plt.figure(figsize=(10, 6))


# Create the primary axis (ax1) for Power Values
ax1 = plt.gca()  # Get the current axis instance

# Plot the timeseries data
#plt.plot(chiller_data['Timestamp'], chiller_data['WKM_eCH45150Leistung_elek_aktuell'], label='Power Data', color='gray')
plt.plot(chiller_data['WKM_eCH45150Leistung_elek_aktuell'], label='Power Data', color='black')
# Plot the power state boundaries
for state, (low, high) in power_states.items():
    #plt.axhline(y=low, color=colors[state], linestyle='--', alpha=0.7, label=f"{state} Lower Bound")
    plt.axhline(y=high, color=colors[state], linestyle='--', alpha=0.7, label=f"{state} Upper Bound")
# Set labels, title, and legend for ax1
#ax1.set_xlabel('Timestamp')
ax1.set_ylabel('Power in kW',fontsize=16)
#ax1.legend(loc='lower right', bbox_to_anchor=(1.05, 1), ncol=1)
#ax1.set_title('Chiller Power Values with Power State Boundaries and EER')

ax1.set_xlabel('Time in s', fontsize=16)
ax1.tick_params(axis='both', which='major', labelsize=16)
# Create a secondary y-axis (ax2) for EER values
ax2 = ax1.twinx()

# Plot the EER values on ax2
ax2.tick_params(axis='both', which='major', labelsize=16)
ax2.plot(chiller_data['WKM_eCH45150EER_aktuell'], label='EER', color='black',linestyle='dotted')
ax2.set_ylabel('EER',fontsize=16)
# Erstellen einer gemeinsamen Legende für beide Achsen unter dem Plot
handles1, labels1 = ax1.get_legend_handles_labels()
handles2, labels2 = ax2.get_legend_handles_labels()
plt.legend(handles1 + handles2, labels1 + labels2, loc='upper center', bbox_to_anchor=(0.5, -0.15), ncol=4, fontsize=12)


plt.xticks(rotation=45)
plt.tight_layout()

#plt.savefig('data_chiller_states.png')
plt.savefig('data_chiller_states.svg',format='svg')
plt.show()

# %%
# CALCULATE HOLDING DURATIONS IN DATASET

import pandas as pd
from datetime import datetime

data_online = pd.read_csv('Merge_chiller_15_310823.csv')

# Definiere die Powerstates und deren Grenzen
power_states = {
    "Powerstate_1": (0, 3),
    "Powerstate_2": (3, 6),
    "Powerstate_3": (6, 9),
    "Powerstate_4": (9, 12),
    "Powerstate_5": (12, 15)
}

# Konvertiere die Timestamps in datetime-Objekte
data_online['Timestamp'] = pd.to_datetime(data_online['Timestamp'])

# Füge neue Spalten für den Powerstate, die Dauer, den Startindex und den Endindex hinzu
data_online['Powerstate'] = None
data_online['Duration'] = 0
data_online['StartIndex'] = None
data_online['EndIndex'] = None

# Initialisiere Variablen für die längste Dauer und die entsprechenden Indizes
longest_duration = {state: {'Duration': 0, 'StartIndex': None, 'EndIndex': None} for state in power_states}

# Bestimme den Powerstate für jeden Messwert und die längste Dauer
current_state = None
for i, row in data_online.iterrows():
    power_value = row['WKM_eCH45150Leistung_elek_aktuell']
    timestamp = row['Timestamp']
    for state_name, (lower_bound, upper_bound) in power_states.items():
        if lower_bound <= power_value < upper_bound:
            data_online.at[i, 'Powerstate'] = state_name
            if state_name == current_state:
                # Fortsetzung des aktuellen Zustands
                duration = (timestamp - data_online.at[i-1, 'Timestamp']).total_seconds()
                data_online.at[i, 'Duration'] = data_online.at[i-1, 'Duration'] + duration
                data_online.at[i, 'StartIndex'] = data_online.at[i-1, 'StartIndex']
                data_online.at[i, 'EndIndex'] = i
            else:
                # Beginn eines neuen Zustands
                data_online.at[i, 'Duration'] = 0
                data_online.at[i, 'StartIndex'] = i
                data_online.at[i, 'EndIndex'] = i
                current_state = state_name
            
            # Aktualisiere die längste Dauer und Indizes falls nötig
            if data_online.at[i, 'Duration'] > longest_duration[state_name]['Duration']:
                longest_duration[state_name]['Duration'] = data_online.at[i, 'Duration']
                longest_duration[state_name]['StartIndex'] = data_online.at[i, 'StartIndex']
                longest_duration[state_name]['EndIndex'] = i
            
            break

# Ausgabe der längsten Dauern und Indizes für jeden Zustand
print("Maximale Dauer für jeden Powerstate in Sekunden und zugehörige Indizes:")
for state, info in longest_duration.items():
    duration = info['Duration']
    start_idx = info['StartIndex']
    end_idx = info['EndIndex']
    print(f"{state}: {duration:.2f} Sekunden, StartIndex: {start_idx}, EndIndex: {end_idx}")


print('>>>>>>>Calculate Holding Durations Successfull<<<<<<')



# %%
# Bestimmung Power gradiente
import numpy as np

power_values = chiller_data['WKM_eCH45150Leistung_elek_aktuell'].values
timestamps = pd.to_datetime(chiller_data['Timestamp'])

# Compute differences between consecutive power values
delta_power = np.diff(power_values)

# Compute differences between consecutive timestamps (in seconds)
delta_time = (timestamps - timestamps.shift()).dropna().dt.total_seconds().values

# Compute rate of change for each interval
rate_of_change = np.abs(delta_power) / delta_time

# Compute mean rate of change
mean_rate_of_change = rate_of_change.mean()

print(f"Mean Rate of Change: {mean_rate_of_change} kW/second")

# %% Calculate maximum of Reaction duration
# max ist von power state 1 zu powert state 5 
T_Reac = (14.5-0.3)/mean_rate_of_change
print(f"Reaction Time max: {T_Reac} second or {T_Reac/60} min")

# Compute rate of change for each interval
rate_of_change = delta_power / delta_time

# Filter out increasing and decreasing rates of change
increasing_rate_of_change = rate_of_change[rate_of_change > 0]
decreasing_rate_of_change = rate_of_change[rate_of_change < 0]

# Calculate max and min for increasing rates of change
max_increasing_roc = increasing_rate_of_change.max()
min_increasing_roc = increasing_rate_of_change.min()

# Calculate max and min for decreasing rates of change
max_decreasing_roc = decreasing_rate_of_change.max()  # This will be closer to zero, as these values are negative
min_decreasing_roc = decreasing_rate_of_change.min()  # This will be the most negative value

print(f"Maximum Rate of Change for Increasing Power Values: {max_increasing_roc} kW/second")
print(f"Minimum Rate of Change for Increasing Power Values: {min_increasing_roc} kW/second")
print(f"Maximum Rate of Change for Decreasing Power Values: {max_decreasing_roc} kW/second")
print(f"Minimum Rate of Change for Decreasing Power Values: {min_decreasing_roc} kW/second")

# %%i
# bestimmung Holding Durations

def calculate_filling_time(energy_capacity_kWh, power_kW):
    """
    Berechnet die Zeit, die benötigt wird, um einen Kältespeicher zu füllen.
    
    Parameter:
    energy_capacity_kWh (float): Die Kapazität des Kältespeichers in kWh.
    power_kW (float): Die Leistung der Kältemaschine in kW.
    
    Rückgabe:
    float: Die Zeit in Sekunden, die benötigt wird, um den Speicher zu füllen.
    """
    if power_kW == 0:
        return float('inf')  # Unendlich, da der Speicher nie gefüllt wird
    else:
        time_seconds = (energy_capacity_kWh / power_kW) * 3600  # Umrechnung in Sekunden
        return time_seconds

# Kapazität des Kältespeichers in kWh
energy_capacity_kWh = 15  

# Leistungszustände in kW (min, max)
power_states = {
    'P1': (0.3, 2.9),
    'P2': (3, 5.9),
    'P3': (6, 8.9),
    'P4': (9, 11.9),
    'P5': (12, 15)
}

# Berechne die minimale und maximale Fülldauer für jeden Leistungszustand
for state, (min_power, max_power) in power_states.items():
    min_time = calculate_filling_time(energy_capacity_kWh, max_power)  # Minimale Fülldauer bei max Leistung
    max_time = calculate_filling_time(energy_capacity_kWh, min_power)  # Maximale Fülldauer bei min Leistung
    print(f"{state}: Minimale Fülldauer = {min_time} Sekunden, Maximale Fülldauer = {max_time} Sekunden")



# %% 
# BUILD final EFDM Set
power_states = {
    "Powerstate_1": (0, 2.9),
    "Powerstate_2": (3, 5.9),
    "Powerstate_3": (6, 8.9),
    "Powerstate_4": (9, 11.9),
    "Powerstate_5": (12, 14.5)
}
print("Maximale Dauer für jeden Powerstate in Sekunden und zugehörige Indizes:")
for state, info in longest_duration.items():
    duration = info['Duration']
    start_idx = info['StartIndex']
    end_idx = info['EndIndex']
    print(f"{state}: {duration:.2f} Sekunden, StartIndex: {start_idx}, EndIndex: {end_idx}")

print(f"Mean Rate of Change: {mean_rate_of_change} kW/second")
print(f"Reaction Time max: {T_Reac} second or {T_Reac/60} min")
print(f"Reaction Time min: {0} second or {0} min")
print("Regenerationszeit konnte über die aufgezeichneten Daten nicht bestimmt werden")



# %% 
#   Analysze Colde Storage Data
m_stor = 1000 # Liter 
cp = 4.184 # 1kJ/(Kg*K)
r_u = 1/3600 # for conversion between J to Wh

T_soll = 15 # °C

T_critic_max = 25 #°C maximale Texhnische Grenze
T_critic_min = 12 #°C minimale technische Grenze
T_flex_max = 20 #°C max Tempe of flex Operation
T_flex_min = T_soll #°C min Temp of flex operation

#Identifizierte Speicherkapazität
import pandas as pd

# Load the dataset
storage_data =  pd.read_csv('Merge_Storage_15_310823.csv')

# Columns to analyze
temperature_columns = ['ThSp_VC100E195ThHy_T_un', 'ThSp_VC100E195ThHy_T_mi', 'ThSp_VC100E195ThHy_T_ob']

# Find min and max values for each column
min_max_values = {}
for column in temperature_columns:
    min_value = storage_data[column].min()
    max_value = storage_data[column].max()
    min_max_values[column] = (min_value, max_value)

print(min_max_values)
import pandas as pd
import matplotlib.pyplot as plt

# Plot the data
plt.figure(figsize=(10, 6))

# for column in temperature_columns:
#     plt.plot(storage_data[column], label=column)

plt.plot(storage_data['ThSp_VC100E195ThHy_T_mi'], label='Storage Temperature', color='black')
plt.xlabel('Time in s',fontsize=16)
plt.ylabel('Temperature in °C',fontsize=16)
ax = plt.gca()  # Greift auf das aktuelle Achsen-Objekt zu
ax.tick_params(axis='both', which='major', labelsize=16)  # Schriftgröße für Haupt-Ticks
#plt.legend()
#plt.title('Temperature Data')
plt.grid(True)
plt.tight_layout()
# plt.savefig('data_cold_storage.png')
plt.savefig('data_cold_storage.svg', format='svg')
plt.show()

# %%
# EFDM Key Figures of Cold Storage
Capacity_storage_max = m_stor*cp*(T_critic_max-T_critic_min)*r_u
print(f"Cooling Storage Capacity: Maximum {Capacity_storage_max} kWh")
Capacity_storage_min = m_stor*cp*(T_critic_min-T_critic_min)*r_u
print(f"Cooling Storage Capacity: Minimum{Capacity_storage_min} kWh")
# Initial Energy Content
E_start = m_stor*cp*(T_flex_max-T_flex_min)*r_u
print(f"Initial Energy Content Maximum {E_start} kWh")

E_target = m_stor*cp*(T_flex_min-T_flex_min)*r_u
print(f"Target Energy Content minimum {E_target} kWh")
E_target = m_stor*cp*(T_flex_max-T_flex_min)*r_u
print(f"Target Energy Content maximum {E_target} kWh")

q_bs = 1.37/24 # kWh/h
# Konvertiere kWh in KJ 
q_bs_neu = 1.37*3600/24 # kJ/h
V_Stor = 1 # m³
roh = 1000 # kg/m³
Delta_T = T_flex_max-T_flex_min # Kelvin
E_loss = q_bs_neu/(V_Stor*cp*roh*Delta_T) # Einhiet 1/h
#Um eine solche Rate in Prozent pro Stunde (%/h) umzurechnen,
# multipliziert man den Wert in [1/h] einfach mit 100, da 1 als 100% ausgedrückt wird.
print(f"Eloss_max = {E_loss*100} %/h") # 

# Calculation of Drain
E_Drain_Factor = 0.3 # Nach Quelle Schraml, Lindner usw
print(f"E_Drain is {E_Drain_Factor}*Power of WZM")

# Supply System 
Supply_system = "Chiller"
print(f"Supplier of storage is {Supply_system}")


# %%   Create Flex Space

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors



def plot_machine(ax, leistung, dauer, aenderungsrate, unabhaengigkeit, cmap, norm):
    heat = plt.cm.gray(norm(unabhaengigkeit))
    #heat = plt.cm.plasma(unabhaengigkeit)
    
    # Plot the top surface
    Z = aenderungsrate * np.ones_like(leistung)
    ax.plot_surface(leistung, dauer, Z, color=heat,shade=False)#, alpha=1)
    
    # Plot the bottom surface
    Z = np.zeros_like(leistung)
    ax.plot_surface(leistung, dauer, Z, color=heat,shade=False, alpha=0.5)
    
    # Plot vertical edges to create a 3D volume
    for i in [0, 1]:
        for j in [0, 1]:
            ax.plot([leistung[i, j], leistung[i, j]], 
                    [dauer[i, j], dauer[i, j]], 
                    [0, aenderungsrate], 
                    color=heat)

# Befülldauer bei mittlerer gemessener Leistung von 5kW des Chillers 
time_hours = 3
# Daten für Chiller + Kältemaschine
leistung = np.array([[0.3, 14.5], [0.3, 14.5]])  # Power State Chiller min max aus handbuch
dauer = np.array([[time_hours, time_hours], [0.5, 0.5]])  # Dauer 
aenderungsrate = 0.005  # Änderungsrate in kW/s
unabhaengigkeit = 0.7  # Grad der Unabhängigkeit

# Erstellen des Plots für Maschine 4
fig = plt.figure(figsize=(10, 10))
cmap = plt.cm.plasma
norm = mcolors.Normalize(vmin=0, vmax=1)

# Plot für Maschine 4
ax = fig.add_subplot(111, projection='3d')
plot_machine(ax, leistung, dauer, aenderungsrate, unabhaengigkeit, cmap, norm)
#ax.set_title("Chiller + Cooling Storage")

# Skalierung anpassen
ax.set_xlim([0, 20]) #Power
ax.set_ylim([0, 5]) # Dauer
ax.set_zlim([0, 0.6]) # Changerate

# Labels setzen
ax.set_xlabel('Power in kW', fontsize=16)
ax.set_ylabel('Duration in h',fontsize=16)
ax.set_zlabel(r'$\nabla$P in kW/s',fontsize=16)
ax.tick_params(axis='both', labelsize=16) 
plt.tight_layout(pad=0.1)

#plt.savefig('flexspace_chiller_storage.png')
plt.savefig('flexspace_chiller_storage.svg', format='svg')
plt.show()
