add version number and remove some draw code

main.py

@@ -1,199 +1,16 @@
 #!/usr/bin/env python
 # coding: utf-8
 
-# In[40]:
+# In[14]:
 
-import pandas as pd
-class pv_config:
-    def __init__(self, capacity, cost_per_kW, lifetime, loss):
-        self.capacity = capacity
-        self.cost_per_kW = cost_per_kW
-        self.lifetime = lifetime
-        self.loss = loss
-    def get_cost(self):
-        return self.capacity * self.cost_per_kW
-    def get_cost_per_year(self):
-        return self.capacity * self.cost_per_kW / self.lifetime
-class ess_config:
-    def __init__(self, capacity, cost_per_kW, lifetime, loss, charge_power, discharge_power):
-        self.capacity = capacity
-        self.cost_per_kW = cost_per_kW
-        self.lifetime = lifetime
-        self.loss = loss
-        self.storage = 100
-        self.charge_power = charge_power
-        self.discharge_power = discharge_power
-    def get_cost(self):
-        return self.capacity * self.cost_per_kW
-    def get_cost_per_year(self):
-        return self.capacity * self.cost_per_kW / self.lifetime
-
-class grid_config:
-    def __init__(self, capacity, grid_loss, sell_price):
-        # self.price_schedule = price_schedule
-        self.loss = grid_loss
-        self.sell_price = sell_price
-        self.capacity = capacity
-    
-    def get_price_for_time(self, time):
-        hour, minute = map(int, time.split(':'))
-        total_minutes = hour * 60 + minute
-        for _, row in self.price_schedule.iterrows():
-            start_hour, start_minute = map(int, row['time_start'].split(':'))
-            end_hour, end_minute = map(int, row['time_end'].split(':'))
-            start_total_minutes = start_hour * 60 + start_minute
-            end_total_minutes = end_hour * 60 + end_minute
-            if start_total_minutes <= total_minutes < end_total_minutes:
-                return row['price']
-        return 0.1  # 默认电价，以防万一没有匹配的时间段
-
-
-
-# In[41]:
-class EnergySystem:
-    def __init__(self, pv_type: pv_config, ess_type: ess_config, grid_type: grid_config):
-        self.pv = pv_type
-        self.ess = ess_type
-        self.grid = grid_type
-        self.day_generated = []
-        self.generated = 0
-        self.stored = 0
-        self.hour_stored = []
-        self.hour_stored_2 = []
-        self.afford = True
-        self.cost = self.ess.get_cost() + self.pv.get_cost()
-        self.overload_cnt = 0
-        self.spring_week_gen = []
-        self.summer_week_gen = []
-        self.autumn_week_gen = []
-        self.winter_week_gen = []
-        self.spring_week_soc = []
-        self.summer_week_soc = []
-        self.autumn_week_soc = []
-        self.winter_week_soc = []
-        self.granularity = 4
-        self.season_step = self.granularity * 24 * 7 * 12
-        self.season_start= self.granularity * 24 * 7 * 2
-        self.week_length = self.granularity * 24 * 7
-        self.unmet = []
-        
-    def get_cost(self):
-        return self.ess.get_cost()+self.pv.get_cost()
-
-    # 优先使用PV供电给工厂 - 如果PV输出能满足工厂的需求，则直接供电，多余的电能用来给ESS充电。
-    # PV不足时使用ESS补充 - 如果PV输出不足以满足工厂需求，首先从ESS获取所需电量。
-    # 如果ESS也不足以满足需求，再从电网获取 - 当ESS中的存储电量也不足以补充时，再从电网购买剩余所需电量。
-    def simulate(self, data, time_interval):
-        total_benefit = 0
-        for index, row in data.iterrows():
-            time = row['time']
-            sunlight_intensity = row['sunlight']
-            factory_demand = row['demand']
-            electricity_price = row['price']
-            # electricity_price = self.grid.get_price_for_time(time)
-
-            if time == '00:00':
-                self.day_generated.append(self.generated)
-                self.generated = 0
-            if time.endswith('14:00'):
-                soc = self.ess.storage / self.ess.capacity
-                self.hour_stored.append(soc)
-            if time.endswith('08:00'):
-                soc = self.ess.storage / self.ess.capacity
-                self.hour_stored_2.append(soc)
-
-            generated_pv_power = self.pv.capacity * sunlight_intensity  # 生成的功率，单位 kW
-            generated_pv_energy = generated_pv_power * time_interval * self.pv.loss  # 生成的能量，单位 kWh
-            self.generated += generated_pv_energy
-            # pv生成的能量如果比工厂的需求要大
-            if generated_pv_energy >= factory_demand * time_interval:
-                # 剩余的能量(kwh) = pv生成的能量 - 工厂需求的功率 * 时间间隔 
-                surplus_energy = generated_pv_energy - factory_demand * time_interval
-                # 要充到ess中的能量 = min(剩余的能量,ess的充电功率*时间间隔(ess在时间间隔内能充进的电量),ess的容量-ess储存的能量(ess中能冲进去的电量))
-                charge_to_ess = min(surplus_energy, self.ess.charge_power * time_interval, self.ess.capacity - self.ess.storage)
-                self.ess.storage += charge_to_ess
-                surplus_after_ess = surplus_energy - charge_to_ess
-                # 如果还有电量盈余,且pv功率大于ess的充电功率+工厂的需求功率则准备卖电
-                if surplus_after_ess > 0 and generated_pv_power > self.ess.charge_power + factory_demand:
-                    sold_to_grid = surplus_after_ess
-                    sell_income = sold_to_grid * self.grid.sell_price
-                    total_benefit += sell_income
-                # 节省的能量 = 工厂需求的能量 * 时间段
-                # total_energy = factory_demand * time_interval
-                saved_energy = factory_demand * time_interval
-            # pv比工厂的需求小
-            else:
-                # 从ess中需要的电量 = 工厂需要的电量 - pv中的电量
-                needed_from_ess = factory_demand * time_interval - generated_pv_energy
-                # 如果ess中存的电量比需要的多
-                if self.ess.storage * self.ess.loss >= needed_from_ess:
-                    # 取出电量
-                    if self.ess.discharge_power * time_interval * self.ess.loss < needed_from_ess:
-                        discharging_power = self.ess.discharge_power * time_interval 
-                    else:
-                        discharging_power = needed_from_ess / self.ess.loss
-
-                    self.ess.storage -= discharging_power
-                    # 节省下来的能量 = pv的能量 + 放出来的能量
-                    saved_energy = generated_pv_energy + discharging_power * self.ess.loss
-                else:
-                    # 如果存的电量不够
-                    # 需要把ess中的所有电量释放出来
-                    if self.grid.capacity * time_interval + generated_pv_energy + self.ess.storage * self.ess.loss < factory_demand * time_interval:
-                        self.afford = False
-                        self.overload_cnt+=1
-                        log = f"index: {index}, time: {time}, SoC:{self.ess.storage / self.ess.capacity}%, storage: {self.ess.storage}, pv_gen:{generated_pv_power}, power_demand: {factory_demand}, overload_cnt:{self.overload_cnt}, day:{int(index/96) + 1}"
-                        self.unmet.append((index,time,factory_demand,generated_pv_power))
-                        # with open(f'plots/summary/ess-{self.ess.capacity}-pv-{self.pv.capacity}', 'a') as f:
-                            # f.write(log)
-                        print(log)
-                        # self.unmet.append(log)
-                    saved_energy = generated_pv_energy + self.ess.storage * self.ess.loss
-                    self.ess.storage = 0
-                    needed_from_grid = factory_demand * time_interval - saved_energy
-                    net_grid = min(self.grid.capacity * time_interval, needed_from_grid) * self.grid.loss
-                    # grid_energy += net_grid
-                    # total_energy += net_grid
-            # print(total_energy)
-            # 工厂需求量-总能量
-            # unmet_demand = max(0, factory_demand * time_interval - total_energy)
-            # benefit = (total_energy - unmet_demand) * electricity_price
-            benefit = (saved_energy) * electricity_price
-            cost = net_grid * electricity_price
-            # print(f"time:{time} benefit: {benefit}, cost: {cost}")
-            total_benefit += benefit - cost
-            # # spring
-            week_start = self.season_start
-            week_end = self.week_length + week_start
-            if index in range(week_start, week_end):
-                self.spring_week_gen.append(generated_pv_power)
-                self.spring_week_soc.append(self.ess.storage / self.ess.capacity)
-            # summer
-            # week_start += self.season_step
-            # week_end += self.season_step
-            # if index in range(week_start, week_end):
-            #     self.summer_week_gen.append(generated_pv_power)
-            #     self.summer_week_soc.append(self.ess.storage / self.ess.capacity)
-            # # autumn
-            # week_start += self.season_step
-            # week_end += self.season_step
-            # if index in range(week_start, week_end):
-            #     self.autumn_week_gen.append(generated_pv_power)
-            #     self.autumn_week_soc.append(self.ess.storage / self.ess.capacity)
-            # week_start += self.season_step
-            # week_end += self.season_step
-            # if index in range(week_start, week_end):
-            #     self.winter_week_gen.append(generated_pv_power)
-            #     self.winter_week_soc.append(self.ess.storage / self.ess.capacity)
-
-        return total_benefit
 
 import os
 import glob
 import shutil
 
 def clear_folder_make_ess_pv(folder_path):
-    shutil.rmtree(folder_path)
+    if os.path.isdir(folder_path):
+        shutil.rmtree(folder_path)
     os.makedirs(folder_path)
     os.makedirs(os.path.join(folder_path,'ess'))
     os.makedirs(os.path.join(folder_path,'pv'))
@@ -202,22 +19,24 @@ folder_path = 'plots'
 clear_folder_make_ess_pv(folder_path)
 
 
-# In[42]:
+# In[15]:
 
 
 import matplotlib.pyplot as plt
 import seaborn as sns
 import numpy as np
 import pandas as pd
+from EnergySystem import EnergySystem
 from config import pv_config, grid_config, ess_config
-figure_size = (10,8)
 
 
-# In[43]:
+# In[16]:
 
 
 import json
 
+print("Version 0.0.2")
+
 with open('config.json', 'r') as f:
     js_data = json.load(f)
 
@@ -243,6 +62,17 @@ pv_groups = js_data["pv_capacities"]["groups"]
 ess_begin = js_data["ess_capacities"]["begin"]
 ess_end = js_data["ess_capacities"]["end"]
 ess_groups = js_data["ess_capacities"]["groups"]
+
+annot_unmet = js_data["annotated"]["unmet_prob"]
+annot_benefit = js_data["annotated"]["benefit"]
+annot_cost = js_data["annotated"]["cost"]
+
+title_unmet = js_data["plot_title"]["unmet_prob"]
+title_cost = js_data["plot_title"]["cost"]
+title_benefit = js_data["plot_title"]["benefit"]
+
+figure_size = (js_data["figure_size"]["length"], js_data["figure_size"]["height"])
+
 pv_capacities = np.linspace(pv_begin, pv_end, pv_groups)
 ess_capacities = np.linspace(ess_begin, ess_end, ess_groups)
 results = pd.DataFrame(index=pv_capacities, columns= ess_capacities)
@@ -251,7 +81,7 @@ costs = pd.DataFrame(index=pv_capacities, columns= ess_capacities)
 overload_cnt = pd.DataFrame(index=pv_capacities, columns= ess_capacities)
 
 
-# In[44]:
+# In[17]:
 
 
 hour_demand = []
@@ -267,7 +97,7 @@ plt.savefig('plots/demand.png')
 plt.close()
 
 
-# In[45]:
+# In[18]:
 
 
 def cal_profit(es: EnergySystem, saved_money):
@@ -275,13 +105,13 @@ def cal_profit(es: EnergySystem, saved_money):
     return profit
 
 
-# In[46]:
+# In[24]:
 
 
-for pv_capacity in pv_capacities:
+for ess_capacity in ess_capacities:
-    print(f"pv_capacity:{pv_capacity}")
+    print(f"ess_capacity:{ess_capacity}")
-    for ess_capacity in ess_capacities:
+    for pv_capacity in pv_capacities:
-        print(f"ess_capacity:{ess_capacity}")
+        print(f"pv_capacity:{ess_capacity}")
         pv = pv_config(capacity=pv_capacity, 
                         cost_per_kW=pv_cost_per_kW,
                         lifetime=pv_lifetime, 
@@ -306,23 +136,23 @@ for pv_capacity in pv_capacities:
         pv_generated = energySystem.day_generated
         ess_generated = energySystem.hour_stored
         ess_generated_2 = energySystem.hour_stored_2
-        plt.figure(figsize=(10,8));
+    plt.figure(figsize=(10,8));
-        plt.plot(ess_generated)
+    plt.plot(ess_generated)
-        plt.xlabel('day #')
+    plt.xlabel('day #')
-        plt.ylabel('SoC %')
+    plt.ylabel('SoC %')
-        plt.title(f'14:00 ESS SoC \n PV cap:{pv_capacity}, ESS cap:{ess_capacity}')
+    plt.title(f'14:00 ESS SoC \n PV cap:{pv_capacity}, ESS cap:{ess_capacity}')
-        plt.savefig(f'plots/ess/1400-{pv_capacity}-{ess_capacity}.png')
+    plt.savefig(f'plots/ess/1400-{pv_capacity}-{ess_capacity}.png')
-        plt.close()
+    plt.close()
-        plt.figure(figsize=(10,8));
+    plt.figure(figsize=(10,8));
-        plt.plot(ess_generated_2)
+    plt.plot(ess_generated_2)
-        plt.xlabel('day #')
+    plt.xlabel('day #')
-        plt.ylabel('SoC%')
+    plt.ylabel('SoC%')
-        plt.title(f'08:00 ESS SoC \n PV cap:{pv_capacity}, ESS cap:{ess_capacity}')
+    plt.title(f'08:00 ESS SoC \n PV cap:{pv_capacity}, ESS cap:{ess_capacity}')
-        plt.savefig(f'plots/ess/0800-{pv_capacity}-{ess_capacity}.png')
+    plt.savefig(f'plots/ess/0800-{pv_capacity}-{ess_capacity}.png')
-        plt.close()
+    plt.close()
-        print(energySystem.unmet)
+        # print(energySystem.unmet)
-        spring_week_start = energySystem.season_start
+        # spring_week_start = energySystem.season_start
-        spring_week_end = spring_week_start + energySystem.week_length
+        # spring_week_end = spring_week_start + energySystem.week_length
         # summer_week_start = energySystem.season_start + 1 * energySystem.season_step
         # summer_week_end = summer_week_start + energySystem.week_length
         # autumn_week_start = energySystem.season_start + 2 * energySystem.season_step
@@ -330,13 +160,13 @@ for pv_capacity in pv_capacities:
         # winter_week_start = energySystem.season_start + 3 * energySystem.season_step
         # winter_week_end = winter_week_start+ energySystem.week_length
 
-        spring_consume_data = []
+        # spring_consume_data = []
         # summer_consume_data = []
         # autumn_consume_data = []
         # winter_consume_data = []
-        for index, row in data.iterrows():
+        # for index, row in data.iterrows():
-            if index in range(spring_week_start, spring_week_end):
+            # if index in range(spring_week_start, spring_week_end):
-                spring_consume_data.append(row['demand'])
+                # spring_consume_data.append(row['demand'])
         #     elif index in range(summer_week_start, summer_week_end):
         #         summer_consume_data.append(row['demand'])
         #     elif index in range(autumn_week_start, autumn_week_end):
@@ -344,12 +174,12 @@ for pv_capacity in pv_capacities:
         #     elif index in range(winter_week_start, winter_week_end):
         #         winter_consume_data.append(row['demand'])
 
-        spring_week_time = list(range(spring_week_start, spring_week_end))
+        # spring_week_time = list(range(spring_week_start, spring_week_end))
         # summer_week_time = list(range(summer_week_start, summer_week_end))
         # autumn_week_time = list(range(autumn_week_start, autumn_week_end))
         # winter_week_time = list(range(winter_week_start, winter_week_end))
 
-        spring_pv_generated = energySystem.spring_week_gen
+        # spring_pv_generated = energySystem.spring_week_gen
         # summer_pv_generated = energySystem.summer_week_gen
         # autumn_pv_generated = energySystem.autumn_week_gen
         # winter_pv_generated = energySystem.winter_week_gen
@@ -362,14 +192,14 @@ for pv_capacity in pv_capacities:
 
         # fig, ax1 = plt.subplots()
 
-        plt.plot(spring_week_time, spring_pv_generated, label = 'pv generation')
+        # plt.plot(spring_week_time, spring_pv_generated, label = 'pv generation')
-        plt.plot(spring_week_time, spring_consume_data, label = 'factory consume')
+        # plt.plot(spring_week_time, spring_consume_data, label = 'factory consume')
-        plt.ylabel('Power / kW')
+        # plt.ylabel('Power / kW')
-        plt.xlabel('15 min #')
+        # plt.xlabel('15 min #')
-        plt.title(f'ess: {energySystem.ess.capacity/1000 } MWh pv: {energySystem.pv.capacity/1000 } MW spring week generate condition')
+        # plt.title(f'ess: {energySystem.ess.capacity/1000 } MWh pv: {energySystem.pv.capacity/1000 } MW spring week generate condition')
-        plt.legend()
+        # plt.legend()
-        plt.savefig(f'plots/{energySystem.ess.capacity}-{energySystem.pv.capacity}-spring.png')
+        # plt.savefig(f'plots/{energySystem.ess.capacity}-{energySystem.pv.capacity}-spring.png')
-        plt.close()
+        # plt.close()
 
         # plt.plot(summer_week_time, summer_pv_generated, label = 'pv generation')
         # plt.plot(summer_week_time, summer_consume_data, label = 'factory consume')
@@ -398,13 +228,13 @@ for pv_capacity in pv_capacities:
         # plt.savefig(f'plots/{energySystem.ess.capacity}-{energySystem.pv.capacity}-winter.png')
         # plt.close()
 
-    plt.figure();
+    # plt.figure();
-    plt.plot(pv_generated)
+    # plt.plot(pv_generated)
-    plt.xlabel('day #')
+    # plt.xlabel('day #')
-    plt.ylabel('Electricity kWh')
+    # plt.ylabel('Electricity kWh')
-    plt.title(f'PV generated pv cap:{pv_capacity}, ess cap:{ess_capacity}')
+    # plt.title(f'PV generated pv cap:{pv_capacity}, ess cap:{ess_capacity}')
-    plt.savefig(f'plots/pv/{pv_capacity}-{ess_capacity}.png')
+    # plt.savefig(f'plots/pv/{pv_capacity}-{ess_capacity}.png')
-    plt.close()
+    # plt.close()
 
 
         # plt.show()
@@ -424,15 +254,25 @@ for pv_capacity in pv_capacities:
     # print(benefit)
 
 
-# In[47]:
+# In[20]:
 
 
-energySystem.unmet
+def save_data(data, filename):
+    data.to_csv(filename+'.csv')
+    data.to_json(filename + '.json')
 
 
-# In[48]:
+# In[21]:
 
 
+import matplotlib.ticker as ticker
+
+if not os.path.isdir('data'):
+    os.makedirs('data')
+
+save_data(results, f'data/{pv_begin}-{pv_end}-{pv_groups}-{ess_begin}-{ess_end}-{ess_groups}-results')
+save_data(costs, f'data/{pv_begin}-{pv_end}-{pv_groups}-{ess_begin}-{ess_end}-{ess_groups}-costs')
+save_data(overload_cnt, f'data/{pv_begin}-{pv_end}-{pv_groups}-{ess_begin}-{ess_end}-{ess_groups}-overload_cnt')
 df=results
 df = df.astype(float)
 df.index = df.index / 1000
@@ -442,15 +282,16 @@ max_value = df.max().max()
 max_scale = max(abs(min_value/1000), abs(max_value/1000))
 plt.figure(figsize=figure_size)
 cmap = sns.color_palette("coolwarm", as_cmap=True)
-sns.heatmap(df/1000, annot=True, fmt=".1f", cmap=cmap, vmin=-max_scale, vmax=max_scale)
+ax = sns.heatmap(df/1000, fmt=".1f", cmap=cmap, vmin=-max_scale, vmax=max_scale, annot=annot_benefit)
-plt.title('Benefit Heatmap Based on PV and ESS Capacities (kEUR/year)')
+# ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%.1f'))
+plt.title(title_benefit)
 plt.gca().invert_yaxis()
 plt.xlabel('ESS Capacity (MWh)')
 plt.ylabel('PV Capacity (MW)')
 plt.savefig('plots/benefit.png')
 
 
-# In[49]:
+# In[22]:
 
 
 df = costs
@@ -459,8 +300,8 @@ df.index = df.index / 1000
 df.columns = df.columns / 1000
 
 plt.figure(figsize=figure_size)
-sns.heatmap(df/1000000, annot=True, fmt=".1f", cmap='viridis')
+sns.heatmap(df/1000000,  fmt=".1f", cmap='viridis', annot=annot_cost)
-plt.title('Costs of the PV System/million Eur')
+plt.title(title_cost)
 plt.gca().invert_yaxis()
 plt.xlabel('ESS Capacity (MWh)')
 plt.ylabel('PV Capacity (MW)')
@@ -475,13 +316,7 @@ plt.savefig('plots/costs.png')
     # print(benefit)
 
 
-# In[ ]:
+# In[23]:
-
-
-
-
-
-# In[50]:
 
 
 from matplotlib.colors import LinearSegmentedColormap
@@ -495,8 +330,13 @@ max_scale = max(abs(min_value/1000), abs(max_value/1000))
 
 plt.figure(figsize=figure_size)
 cmap = LinearSegmentedColormap.from_list("", ["white", "blue"])
-sns.heatmap(df/(4*24*365), annot=True, fmt=".1f", cmap=cmap, vmin=0, vmax=1)
+ax = sns.heatmap(df/(4*24*365), fmt=".00%", cmap=cmap, vmin=0, vmax=1, annot=annot_unmet)
-plt.title('Probability of unmet electricity demands')
+cbar = ax.collections[0].colorbar
+cbar.set_ticks([0, 0.25, 0.5, 0.75, 1])
+cbar.set_ticklabels(['0%', '25%', '50%', '75%', '100%'])
+cbar.ax.yaxis.set_major_formatter(ticker.FuncFormatter(lambda x, pos: f'{x:.0%}'))
+
+plt.title(title_unmet)
 plt.gca().invert_yaxis()
 plt.xlabel('ESS Capacity (MWh)')
 plt.ylabel('PV Capacity (MW)')