#!/usr/bin/env python # coding: utf-8 # In[40]: 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) os.makedirs(folder_path) os.makedirs(os.path.join(folder_path,'ess')) os.makedirs(os.path.join(folder_path,'pv')) folder_path = 'plots' clear_folder_make_ess_pv(folder_path) # In[42]: import matplotlib.pyplot as plt import seaborn as sns import numpy as np import pandas as pd from config import pv_config, grid_config, ess_config figure_size = (10,8) # In[43]: import json with open('config.json', 'r') as f: js_data = json.load(f) data = pd.read_csv('combined_data.csv') time_interval = js_data["time_interval"]["numerator"] / js_data["time_interval"]["denominator"] pv_loss = js_data["pv"]["loss"] pv_cost_per_kW = js_data["pv"]["cost_per_kW"] pv_lifetime = js_data["pv"]["lifetime"] ess_loss = js_data["ess"]["loss"] ess_cost_per_kW = js_data["ess"]["cost_per_kW"] ess_lifetime = js_data["ess"]["lifetime"] grid_loss = js_data["grid"]["loss"] sell_price = js_data["grid"]["sell_price"] #kWh grid_capacity = js_data["grid"]["capacity"] #kWh pv_begin = js_data["pv_capacities"]["begin"] pv_end = js_data["pv_capacities"]["end"] 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"] 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) affords = pd.DataFrame(index=pv_capacities, columns= ess_capacities) costs = pd.DataFrame(index=pv_capacities, columns= ess_capacities) overload_cnt = pd.DataFrame(index=pv_capacities, columns= ess_capacities) # In[44]: hour_demand = [] for index, row in data.iterrows(): time = row['time'] demand = row['demand'] if time.endswith('00'): hour_demand.append(demand) plt.figure(figsize=(10,8)) plt.plot(hour_demand) plt.ylabel('Demand Power / kW') plt.savefig('plots/demand.png') plt.close() # In[45]: def cal_profit(es: EnergySystem, saved_money): profit = saved_money - es.ess.get_cost_per_year() - es.pv.get_cost_per_year() return profit # In[46]: for pv_capacity in pv_capacities: print(f"pv_capacity:{pv_capacity}") for ess_capacity in ess_capacities: print(f"ess_capacity:{ess_capacity}") pv = pv_config(capacity=pv_capacity, cost_per_kW=pv_cost_per_kW, lifetime=pv_lifetime, loss=pv_loss) ess = ess_config(capacity=ess_capacity, cost_per_kW=ess_cost_per_kW, lifetime=ess_lifetime, loss=ess_loss, charge_power=ess_capacity, discharge_power=ess_capacity) grid = grid_config(capacity=grid_capacity, grid_loss=grid_loss, sell_price= sell_price) energySystem = EnergySystem(pv_type=pv, ess_type=ess, grid_type= grid) benefit = energySystem.simulate(data, time_interval) results.loc[pv_capacity,ess_capacity] = cal_profit(energySystem, benefit) affords.loc[pv_capacity,ess_capacity] = energySystem.afford overload_cnt.loc[pv_capacity,ess_capacity] = energySystem.overload_cnt costs.loc[pv_capacity,ess_capacity] = energySystem.ess.capacity * energySystem.ess.cost_per_kW + energySystem.pv.capacity * energySystem.pv.cost_per_kW pv_generated = energySystem.day_generated ess_generated = energySystem.hour_stored ess_generated_2 = energySystem.hour_stored_2 plt.figure(figsize=(10,8)); plt.plot(ess_generated) plt.xlabel('day #') plt.ylabel('SoC %') 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.close() plt.figure(figsize=(10,8)); plt.plot(ess_generated_2) plt.xlabel('day #') plt.ylabel('SoC%') 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.close() print(energySystem.unmet) spring_week_start = energySystem.season_start 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 # autumn_week_end = autumn_week_start + energySystem.week_length # winter_week_start = energySystem.season_start + 3 * energySystem.season_step # winter_week_end = winter_week_start+ energySystem.week_length spring_consume_data = [] # summer_consume_data = [] # autumn_consume_data = [] # winter_consume_data = [] for index, row in data.iterrows(): if index in range(spring_week_start, spring_week_end): 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): # autumn_consume_data.append(row['demand']) # 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)) # 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 # summer_pv_generated = energySystem.summer_week_gen # autumn_pv_generated = energySystem.autumn_week_gen # winter_pv_generated = energySystem.winter_week_gen # spring_soc = energySystem.spring_week_soc # summer_soc = energySystem.summer_week_soc # autumn_soc = energySystem.autumn_week_soc # winter_soc = energySystem.winter_week_soc # fig, ax1 = plt.subplots() plt.plot(spring_week_time, spring_pv_generated, label = 'pv generation') plt.plot(spring_week_time, spring_consume_data, label = 'factory consume') plt.ylabel('Power / kW') plt.xlabel('15 min #') plt.title(f'ess: {energySystem.ess.capacity/1000 } MWh pv: {energySystem.pv.capacity/1000 } MW spring week generate condition') plt.legend() plt.savefig(f'plots/{energySystem.ess.capacity}-{energySystem.pv.capacity}-spring.png') plt.close() # plt.plot(summer_week_time, summer_pv_generated, label = 'pv generation') # plt.plot(summer_week_time, summer_consume_data, label = 'factory consume') # plt.ylabel('Power / kW') # plt.xlabel('15 min #') # plt.title(f'ess: {energySystem.ess.capacity/1000 } MWh pv: {energySystem.pv.capacity/1000 } MW summer week generate condition') # plt.legend() # plt.savefig(f'plots/{energySystem.ess.capacity}-{energySystem.pv.capacity}-summer.png') # plt.close() # plt.plot(autumn_week_time, autumn_pv_generated, label = 'pv generation') # plt.plot(autumn_week_time, autumn_consume_data, label = 'factory consume') # plt.ylabel('Power / kW') # plt.xlabel('15 min #') # plt.title(f'ess: {energySystem.ess.capacity/1000 } MWh pv: {energySystem.pv.capacity/1000 } MW autumn week generate condition') # plt.legend() # plt.savefig(f'plots/{energySystem.ess.capacity}-{energySystem.pv.capacity}-autumn.png') # plt.close() # plt.plot(winter_week_time, winter_pv_generated, label = 'pv generation') # plt.plot(winter_week_time, winter_consume_data, label = 'factory consume') # plt.ylabel('Power / kW') # plt.xlabel('15 min #') # plt.title(f'ess: {energySystem.ess.capacity/1000 } MWh pv: {energySystem.pv.capacity/1000 } MW winter week generate condition') # plt.legend() # plt.savefig(f'plots/{energySystem.ess.capacity}-{energySystem.pv.capacity}-winter.png') # plt.close() plt.figure(); plt.plot(pv_generated) plt.xlabel('day #') plt.ylabel('Electricity kWh') plt.title(f'PV generated pv cap:{pv_capacity}, ess cap:{ess_capacity}') plt.savefig(f'plots/pv/{pv_capacity}-{ess_capacity}.png') plt.close() # plt.show() # results = results.astype(float) # pv = pv_config(capacity=100000,cost_per_kW=200,lifetime=25,loss=0.95) # ess = ess_config(capacity=100000,cost_per_kW=300,lifetime=25,loss=0.95,charge_power=100000,discharge_power=100000) # grid = grid_config(price_schedule=price_schedule, capacity=5000, grid_loss=0.95, sell_price=0.4) # grid = grid_config(capacity=50000, grid_loss=0.95, sell_price=0.4) # print(benefit) # In[47]: energySystem.unmet # In[48]: df=results df = df.astype(float) df.index = df.index / 1000 df.columns = df.columns / 1000 min_value = df.min().min() 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) plt.title('Benefit Heatmap Based on PV and ESS Capacities (kEUR/year)') plt.gca().invert_yaxis() plt.xlabel('ESS Capacity (MWh)') plt.ylabel('PV Capacity (MW)') plt.savefig('plots/benefit.png') # In[49]: df = costs df = df.astype(int) 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') plt.title('Costs of the PV System/million Eur') plt.gca().invert_yaxis() plt.xlabel('ESS Capacity (MWh)') plt.ylabel('PV Capacity (MW)') plt.savefig('plots/costs.png') # pv = pv_config(capacity=100000,cost_per_kW=200,lifetime=25,loss=0.95) # ess = ess_config(capacity=100000,cost_per_kW=300,lifetime=25,loss=0.95,charge_power=100000,discharge_power=100000) # grid = grid_config(price_schedule=price_schedule, capacity=5000, grid_loss=0.95, sell_price=0.4) # grid = grid_config(capacity=50000, grid_loss=0.95, sell_price=0.4) # print(benefit) # In[ ]: # In[50]: from matplotlib.colors import LinearSegmentedColormap df = overload_cnt df = df.astype(int) df.index = df.index / 1000 df.columns = df.columns / 1000 min_value = df.min().min() max_value = df.max().max() 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) plt.title('Probability of unmet electricity demands') plt.gca().invert_yaxis() plt.xlabel('ESS Capacity (MWh)') plt.ylabel('PV Capacity (MW)') plt.savefig('plots/unmet.png')