Skip to content

Commit

Permalink
Merge pull request #7 from pjstanle/h2_dispatch
Browse files Browse the repository at this point in the history 
PEM electrolyzer and dispatch
  • Loading branch information
@barker59
barker59 authored Sep 21, 2021
2 parents f3f9024 + f57ddc6 commit c06c8a3
Show file tree
Hide file tree
Showing 6 changed files with 282 additions and 116 deletions.
163 changes: 117 additions & 46 deletions examples/H2 Analysis/h2_analysis_H2A_Refactor.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
from hybrid.solar_source import SolarPlant
from hybrid.wind_source import WindPlant
import PySAM.Singleowner as so
from matplotlib import use
import pandas as pd
import pickle
import json
Expand All @@ -19,8 +20,16 @@
from hopp_for_h2 import hopp_for_h2
from run_h2a import run_h2a
from simple_dispatch import SimpleDispatch
import run_h2_PEM
import numpy as np
from lcoe.lcoe import lcoe as lcoe_calc
import matplotlib.pyplot as plt




import warnings
warnings.filterwarnings("ignore")

# Set API key
load_dotenv()
Expand Down Expand Up @@ -114,14 +123,13 @@ def establish_save_output_dict():
storage_used = False
battery_can_grid_charge = False
grid_connected_hopp = True
kw_continuous = 5000
electrolyzer_sizes = [50, 100, 150, 200]
interconnection_size_mw = 150
load = [kw_continuous for x in
range(0, 8760)] # * (sin(x) + pi) Set desired/required load profile for plant
electrolyzer_sizes = [50, 100, 150, 200]


# Step 2: Load scenarios from .csv and enumerate
scenarios_df = pd.read_csv('H2 Baseline Future Scenarios Test Refactor.csv')
# scenarios_df = pd.read_csv('H2 Baseline Future Scenarios Test Refactor.csv')
scenarios_df = pd.read_csv('single_scenario.csv')
for electrolyzer_size in electrolyzer_sizes:
for critical_load_factor in critical_load_factor_list:
for i, scenario in scenarios_df.iterrows():
Expand All @@ -131,7 +139,12 @@ def establish_save_output_dict():
# -Pass through rotor diameter to pySAM
# -Add wind, solar, storage installed costs
# -Fix "H2 H2 xxx" text
print(scenario)
# print(scenario)

kw_continuous = electrolyzer_size*1000
load = [kw_continuous for x in
range(0, 8760)] # * (sin(x) + pi) Set desired/required load profile for plant

scenario_choice = scenario['Scenario Name']
site_name = scenario['Site Name']
sample_site['lat'] = scenario['Lat']
Expand Down Expand Up @@ -161,6 +174,10 @@ def establish_save_output_dict():
storage_cost_kwh = scenario['Storage Cost kWh']
debt_equity_split = scenario['Debt Equity']

buy_price = scenario['Buy From Grid ($/kWh)']
sell_price = scenario['Sell To Grid ($/kWh)']


#Todo: Add useful life to .csv scenario input instead
scenario['Useful Life'] = useful_life

Expand Down Expand Up @@ -207,7 +224,38 @@ def establish_save_output_dict():
storage_size_mw * 1000)

battery_used, excess_energy, battery_SOC = bat_model.run()
combined_pv_wind_storage_power_production_hopp = combined_pv_wind_power_production_hopp + excess_energy
combined_pv_wind_storage_power_production_hopp = combined_pv_wind_power_production_hopp + battery_used

if sell_price:
profit_from_selling_to_grid = np.sum(excess_energy)*sell_price
else:
profit_from_selling_to_grid = 0.0

if buy_price:
cost_to_buy_from_grid = 0.0

for i in range(len(combined_pv_wind_storage_power_production_hopp)):
if combined_pv_wind_storage_power_production_hopp[i] < kw_continuous:
cost_to_buy_from_grid += (kw_continuous-combined_pv_wind_storage_power_production_hopp[i])*buy_price
combined_pv_wind_storage_power_production_hopp[i] = kw_continuous
else:
cost_to_buy_from_grid = 0.0

plot_battery = False
plot_power = True

if plot_battery:
plt.plot(excess_energy,color="C1",label="excess_energy")
plt.plot(battery_SOC,"--",color="C2",label="battery_SOC")
plt.plot(battery_used,color="C0",label="battery_used")
plt.legend()
plt.show()

if plot_power:
plt.plot(combined_pv_wind_power_production_hopp,label="without battery")
plt.plot(combined_pv_wind_storage_power_production_hopp,"--",label="with battery")
plt.legend()
plt.show()

# Step 6: Run the Python H2A model
# ------------------------- #
Expand All @@ -217,10 +265,21 @@ def establish_save_output_dict():

electrical_generation_timeseries = combined_pv_wind_storage_power_production_hopp

H2_Results, H2A_Results = run_h2a(electrical_generation_timeseries, kw_continuous, electrolyzer_size,
hybrid_plant, reopt_results, scenario,
combined_pv_wind_curtailment_hopp, lcoe, force_electrolyzer_cost, forced_electrolyzer_cost,
total_system_electrical_usage=55.5)
# Old way
# H2_Results, H2A_Results = run_h2a(electrical_generation_timeseries, kw_continuous, electrolyzer_size,
# hybrid_plant, reopt_results, scenario,
# combined_pv_wind_curtailment_hopp, lcoe, force_electrolyzer_cost, forced_electrolyzer_cost,
# total_system_electrical_usage=55.5)

# Parangat model
adjusted_installed_cost = hybrid_plant.grid.financial_model.Outputs.adjusted_installed_cost
useful_life = scenario['Useful Life']
net_capital_costs = reopt_results['outputs']['Scenario']['Site'] \
['Financial']['net_capital_costs']

H2_Results, H2A_Results = run_h2_PEM.run_h2_PEM(electrical_generation_timeseries,turbine_rating,electrolyzer_size,
kw_continuous,forced_electrolyzer_cost,lcoe,adjusted_installed_cost,useful_life,
net_capital_costs)

# Step 6.5: Intermediate financial calculation
#TODO:
Expand All @@ -243,48 +302,57 @@ def establish_save_output_dict():
total_system_installed_cost = total_hopp_installed_cost + total_electrolyzer_cost
annual_operating_cost_hopp = (wind_size_mw * 1000 * 42) + (solar_size_mw * 1000 * 13)
annual_operating_cost_h2 = H2A_Results['Fixed O&M'] * H2_Results['hydrogen_annual_output']
total_annual_operating_costs = annual_operating_cost_hopp + annual_operating_cost_h2
total_annual_operating_costs = annual_operating_cost_hopp + annual_operating_cost_h2 + cost_to_buy_from_grid - profit_from_selling_to_grid
h_lcoe_no_op_cost = lcoe_calc((H2_Results['hydrogen_annual_output']), total_system_installed_cost,
0, 0.07, useful_life)

h_lcoe = lcoe_calc((H2_Results['hydrogen_annual_output']), total_system_installed_cost,
total_annual_operating_costs, 0.07, useful_life)

# Step 7: Print results
# ------------------------- #
#TODO: Tidy up these print statements
print("Future Scenario: {}".format(scenario['Scenario Name']))
print("Wind Cost per KW: {}".format(scenario['Wind Cost KW']))
print("PV Cost per KW: {}".format(scenario['Solar Cost KW']))
print("Storage Cost per KW: {}".format(scenario['Storage Cost kW']))
print("Storage Cost per KWh: {}".format(scenario['Storage Cost kWh']))
print("Wind Size built: {}".format(wind_size_mw))
print("PV Size built: {}".format(solar_size_mw))
print("Storage Size built: {}".format(storage_size_mw))
print("Storage Size built: {}".format(storage_size_mwh))
print("Levelized cost of Electricity (HOPP): {}".format(lcoe))
print("Total Yearly Electrical Output: {}".format(total_elec_production))
print("Total Yearly Hydrogen Production: {}".format(H2_Results['hydrogen_annual_output']))
print("Levelized Cost H2/kg (new method - no operational costs)".format(h_lcoe_no_op_cost))
print("Capacity Factor of Electrolyzer: {}".format(H2_Results['cap_factor']))

print("Levelized cost of H2 (electricity feedstock) (HOPP): {}".format(
H2_Results['feedstock_cost_h2_levelized_hopp']))
print("Levelized cost of H2 (excl. electricity) (H2A): {}".format(H2A_Results['Total Hydrogen Cost ($/kgH2)']))
print("Total unit cost of H2 ($/kg) : {}".format(H2_Results['total_unit_cost_of_hydrogen']))
print("kg H2 cost from net cap cost/lifetime h2 production (HOPP): {}".format(
H2_Results['feedstock_cost_h2_via_net_cap_cost_lifetime_h2_hopp']))

print_reults = False
print_h2_results = True

if print_reults:
# ------------------------- #
#TODO: Tidy up these print statements
print("Future Scenario: {}".format(scenario['Scenario Name']))
print("Wind Cost per KW: {}".format(scenario['Wind Cost KW']))
print("PV Cost per KW: {}".format(scenario['Solar Cost KW']))
print("Storage Cost per KW: {}".format(scenario['Storage Cost kW']))
print("Storage Cost per KWh: {}".format(scenario['Storage Cost kWh']))
print("Wind Size built: {}".format(wind_size_mw))
print("PV Size built: {}".format(solar_size_mw))
print("Storage Size built: {}".format(storage_size_mw))
print("Storage Size built: {}".format(storage_size_mwh))
print("Levelized cost of Electricity (HOPP): {}".format(lcoe))
print("Total Yearly Electrical Output: {}".format(total_elec_production))
print("Total Yearly Hydrogen Production: {}".format(H2_Results['hydrogen_annual_output']))
print("Levelized Cost H2/kg (new method - no operational costs)".format(h_lcoe_no_op_cost))
print("Capacity Factor of Electrolyzer: {}".format(H2_Results['cap_factor']))

if print_h2_results:
# print("Levelized cost of H2 (electricity feedstock) (HOPP): {}".format(
# H2_Results['feedstock_cost_h2_levelized_hopp']))
# print("Levelized cost of H2 (excl. electricity) (H2A): {}".format(H2A_Results['Total Hydrogen Cost ($/kgH2)']))
# print("Total unit cost of H2 ($/kg) : {}".format(H2_Results['total_unit_cost_of_hydrogen']))
# print("kg H2 cost from net cap cost/lifetime h2 production (HOPP): {}".format(
# H2_Results['feedstock_cost_h2_via_net_cap_cost_lifetime_h2_hopp']))
print("h_lcoe: ", h_lcoe)

# Step 8: Plot REopt results
plot_reopt_results(REoptResultsDF, site_name, atb_year, critical_load_factor,
useful_life, tower_height,
wind_size_mw, solar_size_mw, storage_size_mw, storage_size_mwh, lcoe,
H2_Results['feedstock_cost_h2_via_net_cap_cost_lifetime_h2_hopp'],
H2_Results['feedstock_cost_h2_levelized_hopp'],
hybrid_installed_cost, H2A_Results['Total Hydrogen Cost ($/kgH2)'],
H2_Results['total_unit_cost_of_hydrogen'],
output_dir='results/',
monthly_separation=False, reopt_was_run=run_reopt_flag)
plot_reopt = False
if plot_reopt:
plot_reopt_results(REoptResultsDF, site_name, atb_year, critical_load_factor,
useful_life, tower_height,
wind_size_mw, solar_size_mw, storage_size_mw, storage_size_mwh, lcoe,
H2_Results['feedstock_cost_h2_via_net_cap_cost_lifetime_h2_hopp'],
H2_Results['feedstock_cost_h2_levelized_hopp'],
hybrid_installed_cost, H2A_Results['Total Hydrogen Cost ($/kgH2)'],
H2_Results['total_unit_cost_of_hydrogen'],
output_dir='results/',
monthly_separation=False, reopt_was_run=run_reopt_flag)

# Step 9: Plot HOPP Production, Curtailment, and Hydrogen Production Profiles

Expand Down Expand Up @@ -345,5 +413,8 @@ def establish_save_output_dict():
# save_all_runs = save_all_runs.append(save_outputs_dict, sort=False)

# Create dataframe from outputs and save
save_outputs_dict_df = pd.DataFrame(save_outputs_dict)
save_outputs_dict_df.to_csv("results/H2_Analysis_{}.csv".format(site_name))

save_outputs = False
if save_outputs:
save_outputs_dict_df = pd.DataFrame(save_outputs_dict)
save_outputs_dict_df.to_csv("results/H2_Analysis_{}_2.csv".format(site_name))
7 changes: 4 additions & 3 deletions examples/H2 Analysis/hopp_for_h2.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,14 +75,15 @@ def hopp_for_h2(site, scenario, technologies, wind_size_mw, solar_size_mw, stora
energy_shortfall_hopp = [y - x for x, y in
zip(hybrid_plant.grid.generation_profile_from_system[0:8759], load)]
energy_shortfall_hopp = [x if x > 0 else 0 for x in energy_shortfall_hopp]
combined_pv_wind_power_production_hopp = hybrid_plant.grid.system_model.Outputs.system_pre_interconnect_kwac[0:8759]
# combined_pv_wind_power_production_hopp = hybrid_plant.grid.system_model.Outputs.system_pre_interconnect_kwac[0:8759]
combined_pv_wind_power_production_hopp = hybrid_plant.grid.system_model.Outputs.gen[0:8759]
combined_pv_wind_curtailment_hopp = [x - y for x, y in zip(
hybrid_plant.grid.system_model.Outputs.system_pre_interconnect_kwac[0:8759],
hybrid_plant.grid.system_model.Outputs.gen[0:8759])]

# super simple dispatch battery model with no forecasting TODO: add forecasting
print("Length of 'energy_shortfall_hopp is {}".format(len(energy_shortfall_hopp)))
print("Length of 'combined_pv_wind_curtailment_hopp is {}".format(len(combined_pv_wind_curtailment_hopp)))
# print("Length of 'energy_shortfall_hopp is {}".format(len(energy_shortfall_hopp)))
# print("Length of 'combined_pv_wind_curtailment_hopp is {}".format(len(combined_pv_wind_curtailment_hopp)))
# TODO: Fix bug in dispatch model that errors when first curtailment >0
combined_pv_wind_curtailment_hopp[0] = 0

Expand Down
88 changes: 88 additions & 0 deletions examples/H2 Analysis/run_h2_PEM.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,88 @@
from hybrid.PEM_H2_LT_electrolyzer import PEM_electrolyzer_LT
from numpy.lib.function_base import average
import H2AModel
import numpy as np
import pandas as pd


def run_h2_PEM(electrical_generation_timeseries, turbine_rating, electrolyzer_size,
kw_continuous,forced_electrolyzer_cost_kw,lcoe,
adjusted_installed_cost,useful_life,net_capital_costs,
voltage_type="constant", stack_input_voltage_DC=250, min_V_cell=1.62,
p_s_h2_bar=31, stack_input_current_lower_bound=500, cell_active_area=1250,
N_cells=130, total_system_electrical_usage=55.5):

in_dict = dict()
out_dict = dict()
el = PEM_electrolyzer_LT(in_dict, out_dict,electrical_generation_timeseries)

el.power_supply_rating_MW = turbine_rating
el.input_dict['voltage_type'] = voltage_type
el.stack_input_voltage_DC = stack_input_voltage_DC
# Assumptions:
el.min_V_cell = min_V_cell # Only used in variable voltage scenario
el.p_s_h2_bar = p_s_h2_bar # H2 outlet pressure
el.stack_input_current_lower_bound = stack_input_current_lower_bound
el.stack_rating_kW = electrolyzer_size # 1 MW
el.cell_active_area = cell_active_area
el.N_cells = N_cells


el.h2_production_rate()


avg_generation = np.mean(electrical_generation_timeseries) # Avg Generation
# print("avg_generation: ", avg_generation)
cap_factor = avg_generation / kw_continuous

hydrogen_hourly_production = out_dict['h2_produced_kg_hr_system']
print("cap_factor: ", cap_factor)

# Get Daily Hydrogen Production - Add Every 24 hours
i = 0
daily_H2_production = []
while i < 8760:
x = sum(hydrogen_hourly_production[i:i + 25])
daily_H2_production.append(x)
i = i + 25

avg_daily_H2_production = np.mean(daily_H2_production) # kgH2/day
hydrogen_annual_output = sum(hydrogen_hourly_production) # kgH2/year
# elec_remainder_after_h2 = combined_pv_wind_curtailment_hopp

H2A_Results = H2AModel.H2AModel(cap_factor, avg_daily_H2_production, hydrogen_annual_output, force_system_size=True,
forced_system_size=electrolyzer_size, force_electrolyzer_cost=True,
forced_electrolyzer_cost_kw=forced_electrolyzer_cost_kw)


feedstock_cost_h2_levelized_hopp = lcoe * total_system_electrical_usage / 100 # $/kg
# Hybrid Plant - levelized H2 Cost - HOPP
feedstock_cost_h2_via_net_cap_cost_lifetime_h2_hopp = adjusted_installed_cost / \
(hydrogen_annual_output * useful_life) # $/kgH2

# Total Hydrogen Cost ($/kgH2)
h2a_costs = H2A_Results['Total Hydrogen Cost ($/kgH2)']
total_unit_cost_of_hydrogen = h2a_costs + feedstock_cost_h2_levelized_hopp
feedstock_cost_h2_via_net_cap_cost_lifetime_h2_reopt = net_capital_costs / (
(kw_continuous / total_system_electrical_usage) * (8760 * useful_life))

H2_Results = {'hydrogen_annual_output':
hydrogen_annual_output,
'feedstock_cost_h2_levelized_hopp':
feedstock_cost_h2_levelized_hopp,
'feedstock_cost_h2_via_net_cap_cost_lifetime_h2_hopp':
feedstock_cost_h2_via_net_cap_cost_lifetime_h2_hopp,
'feedstock_cost_h2_via_net_cap_cost_lifetime_h2_reopt':
feedstock_cost_h2_via_net_cap_cost_lifetime_h2_reopt,
'total_unit_cost_of_hydrogen':
total_unit_cost_of_hydrogen,
'cap_factor':
cap_factor
}

return H2_Results, H2A_Results





Loading

0 comments on commit c06c8a3

Please sign in to comment.