Merge pull request #8 from pjstanle/develop

Develop

examples/H2 Analysis/H2AModel.py

@@ -853,7 +853,7 @@ def npv(rate, values):
                                      'LCOH Cost Contribution', 'Tax Incentives']) / final_data.loc[
                                     'LCOH Cost Contribution', 'H2 Sales (kg)']) * (
                                        1 + inflation_rate) ** length_of_construction_period / inflation_factor
-    print(Final_Hydrogen_Cost_Real)
+    # print(Final_Hydrogen_Cost_Real)
 
     Cost_Breakdown = pd.DataFrame(columns=['After Tax Present Value', '% of Total', '$/kg of H2'])
     Cost_Breakdown.loc['Capital Related Costs', 'After Tax Present Value'] = -(((df.loc[
@@ -894,7 +894,7 @@ def npv(rate, values):
     # final_data.to_csv('LCOH_Cost_Contribution')
     # Cost_Breakdown.to_csv('Cost_Breakdown.csv')
 
-    print(Cost_Breakdown)
+    # print(Cost_Breakdown)
     feedstock_cost_h2_levelized = Cost_Breakdown.loc['Other Variable Costs (Utilities)', '$/kg of H2']
     results = dict()
     results['Capital Related Costs'] = Cost_Breakdown.loc['Capital Related Costs', '$/kg of H2']

examples/H2 Analysis/h2_analysis_H2A_Refactor.py

@@ -126,6 +126,13 @@ def establish_save_output_dict():
 interconnection_size_mw = 150
 electrolyzer_sizes = [50, 100, 150, 200]
 
+# which plots to show
+plot_power_production = False
+plot_battery = False
+plot_grid = False
+plot_h2 = True
+plot_reopt = False
+
 
 # Step 2: Load scenarios from .csv and enumerate
 # scenarios_df = pd.read_csv('H2 Baseline Future Scenarios Test Refactor.csv')
@@ -218,19 +225,55 @@ def establish_save_output_dict():
             solar_installed_cost = hybrid_plant.solar.financial_model.SystemCosts.total_installed_cost
             hybrid_installed_cost = hybrid_plant.grid.financial_model.SystemCosts.total_installed_cost
 
+            if plot_power_production:
+                plt.title("HOPP power production")
+                plt.plot(combined_pv_wind_power_production_hopp[0:100],label="wind + pv")
+                plt.plot(energy_shortfall_hopp[0:100],label="shortfall")
+                plt.plot(combined_pv_wind_curtailment_hopp[0:100],label="curtailment")
+                plt.plot(load[0:100],label="electrolyzer rating")
+                plt.legend()
+                plt.show()
+
             # Step 5: Run Simple Dispatch Model
             # ------------------------- #
-            bat_model = SimpleDispatch(combined_pv_wind_curtailment_hopp, energy_shortfall_hopp, len(energy_shortfall_hopp),
-                                       storage_size_mw * 1000)
+
+            bat_model = SimpleDispatch()
+            bat_model.Nt = len(energy_shortfall_hopp)
+            bat_model.curtailment = combined_pv_wind_curtailment_hopp
+            bat_model.shortfall = energy_shortfall_hopp
+            bat_model.size_battery = 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 + battery_used
 
+            if plot_battery:
+                plt.figure(figsize=(6,3))
+                plt.subplot(121)
+                plt.plot(combined_pv_wind_curtailment_hopp[0:100],label="curtailment")
+                plt.plot(energy_shortfall_hopp[0:100],label="shortfall")
+                plt.plot(battery_SOC[0:100],label="state of charge")
+                # plt.plot(excess_energy[0:100],label="excess")
+                plt.plot(battery_used[0:100],"--",label="battery used")
+                plt.legend()
+
+                plt.subplot(122)
+                plt.plot(combined_pv_wind_storage_power_production_hopp[0:100],label="wind+pv+storage")
+                plt.plot(combined_pv_wind_power_production_hopp[0:100],"--",label="wind+pv")
+                plt.plot(load[0:100],"--",label="electrolyzer rating")
+
+                plt.legend()
+                plt.tight_layout()
+                plt.show()
+
+            if plot_grid:
+                plt.plot(combined_pv_wind_storage_power_production_hopp[0:100],label="before buy from grid")
+
             if sell_price:
                 profit_from_selling_to_grid = np.sum(excess_energy)*sell_price
             else:
                 profit_from_selling_to_grid = 0.0
 
+            buy_price = False # if you want to force no buy from grid
             if buy_price:
                 cost_to_buy_from_grid = 0.0
 
@@ -241,29 +284,25 @@ def establish_save_output_dict():
             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")
+            energy_to_electrolyzer = [x if x < kw_continuous else kw_continuous for x in combined_pv_wind_storage_power_production_hopp]
+
+            if plot_grid:
+                plt.plot(combined_pv_wind_storage_power_production_hopp[0:100],"--",label="after buy from grid")
+                plt.plot(energy_to_electrolyzer[0:100],"--",label="energy to electrolyzer")
                 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
             # ------------------------- #
             #TODO: Refactor H2A model call
             # Should take as input (electrolyzer size, cost, electrical timeseries, total system electrical usage (kwh/kg),
             # Should give as ouptut (h2 costs by net cap cost, levelized, total_unit_cost of hydrogen etc)   )
 
-            electrical_generation_timeseries = combined_pv_wind_storage_power_production_hopp
+
+            # electrical_generation_timeseries = combined_pv_wind_storage_power_production_hopp
+            electrical_generation_timeseries = np.zeros_like(energy_to_electrolyzer)
+            electrical_generation_timeseries[:] = energy_to_electrolyzer[:]
 
             # Old way
             # H2_Results, H2A_Results = run_h2a(electrical_generation_timeseries, kw_continuous, electrolyzer_size,
@@ -277,10 +316,36 @@ def establish_save_output_dict():
             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,
+            # intalled costs:
+            # hybrid_plant.grid.financial_model.costs
+
+            # system_rating = electrolyzer_size
+            system_rating = wind_size_mw + solar_size_mw
+            H2_Results, H2A_Results = run_h2_PEM.run_h2_PEM(electrical_generation_timeseries,electrolyzer_size,
                             kw_continuous,forced_electrolyzer_cost,lcoe,adjusted_installed_cost,useful_life,
                             net_capital_costs)
 
+            if plot_h2:
+                hydrogen_hourly_production = H2_Results['hydrogen_hourly_production']
+                plt.figure(figsize=(6,3))
+
+                plt.subplot(121)
+                plt.plot(electrical_generation_timeseries[0:100])
+                plt.ylim(0,max(electrical_generation_timeseries[0:100])*1.2)
+                plt.plot(load[0:100],label="electrolyzer rating")
+                plt.title("energy to electrolyzer")
+
+                plt.subplot(122)
+                plt.plot(hydrogen_hourly_production[0:100])
+                plt.ylim(0,max(hydrogen_hourly_production[0:100])*1.2)
+                plt.title("hydrogen production")
+
+                plt.tight_layout()
+                plt.show()
+
+
+
+
             # Step 6.5: Intermediate financial calculation
             #TODO:
             # - Get Hybrid installed cost (wind, solar, storage)
@@ -342,7 +407,7 @@ def establish_save_output_dict():
                 print("h_lcoe: ", h_lcoe)
 
             # Step 8: Plot REopt results
-            plot_reopt = False
+
             if plot_reopt:
                 plot_reopt_results(REoptResultsDF, site_name, atb_year, critical_load_factor,
                                 useful_life, tower_height,

examples/H2 Analysis/hopp_for_h2.py

@@ -72,14 +72,13 @@ def hopp_for_h2(site, scenario, technologies, wind_size_mw, solar_size_mw, stora
     hybrid_plant.simulate(scenario['Useful Life'])
 
     # HOPP Specific Energy Metrics
-    energy_shortfall_hopp = [y - x for x, y in
-                             zip(hybrid_plant.grid.generation_profile_from_system[0:8759], load)]
+    combined_pv_wind_power_production_hopp = hybrid_plant.grid.system_model.Outputs.system_pre_interconnect_kwac[0:8759]
+    energy_shortfall_hopp = [x - y for x, y in
+                             zip(load,combined_pv_wind_power_production_hopp)]
     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.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])]
+    combined_pv_wind_curtailment_hopp = [x - y for x, y in
+                             zip(combined_pv_wind_power_production_hopp,load)]
+    combined_pv_wind_curtailment_hopp = [x if x > 0 else 0 for x in combined_pv_wind_curtailment_hopp]
 
     # super simple dispatch battery model with no forecasting TODO: add forecasting
     # print("Length of 'energy_shortfall_hopp is {}".format(len(energy_shortfall_hopp)))

examples/H2 Analysis/run_h2_PEM.py

@@ -5,7 +5,7 @@
 import pandas as pd
 
 
-def run_h2_PEM(electrical_generation_timeseries, turbine_rating, electrolyzer_size,
+def run_h2_PEM(electrical_generation_timeseries, 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,
@@ -16,18 +16,22 @@ def run_h2_PEM(electrical_generation_timeseries, turbine_rating, electrolyzer_si
     out_dict = dict()
     el = PEM_electrolyzer_LT(in_dict, out_dict,electrical_generation_timeseries)
 
-    el.power_supply_rating_MW = turbine_rating
+    # el.power_supply_rating_MW = electrolyzer_size
+    # el.power_supply_rating_MW = power_supply_rating_MW
+    print("electrolyzer size: ", electrolyzer_size)
+    el.electrolyzer_system_size_MW = electrolyzer_size
     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
 
 
+    print("running production rate")
     el.h2_production_rate()
 
 
@@ -36,7 +40,7 @@ def run_h2_PEM(electrical_generation_timeseries, turbine_rating, electrolyzer_si
     cap_factor = avg_generation / kw_continuous
 
     hydrogen_hourly_production = out_dict['h2_produced_kg_hr_system']
-    print("cap_factor: ", cap_factor)
+    # print("cap_factor: ", cap_factor)
 
     # Get Daily Hydrogen Production - Add Every 24 hours
     i = 0
@@ -70,14 +74,16 @@ def run_h2_PEM(electrical_generation_timeseries, turbine_rating, electrolyzer_si
                         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_hopp,
-                   'feedstock_cost_h2_via_net_cap_cost_lifetime_h2_reopt':
+                    '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':
                        total_unit_cost_of_hydrogen,
-                   'cap_factor':
-                       cap_factor
+                    'cap_factor':
+                       cap_factor,
+                    'hydrogen_hourly_production':
+                        hydrogen_hourly_production
                    }
 
     return H2_Results, H2A_Results

examples/H2 Analysis/simple_dispatch.py

@@ -4,23 +4,20 @@
 
 class SimpleDispatch():
 
-    def __init__(self, combined_pv_wind_curtailment_hopp, energy_shortfall_hopp, N, size_battery):
+    def __init__(self):
 
+        # length of simulation
+        self.Nt = 1
+
         # amount of curtailment experienced by plant
-        self.curtailment = combined_pv_wind_curtailment_hopp
+        self.curtailment = np.zeros(self.Nt)
 
         # amount of energy needed from the battery
-        self.shortfall = energy_shortfall_hopp
-        # print("Energy shortfall in battery simulation is: {}".format(energy_shortfall_hopp))
-
-        # length of simulation
-        self.Nt = N
-        # print("Length of battery simulation is: {}".format(N))
-
+        self.shortfall = np.zeros(self.Nt)
+
         # size of battery (assumed to be the same as the charge rate per hour)
-        self.size_battery = size_battery
-        # print("Size Battery in battery simulation is: {}".format(size_battery))
-
+        self.size_battery = 0
+
 
     def run(self):
 

hybrid/PEM_H2_LT_electrolyzer.py

@@ -60,7 +60,10 @@ def __init__(self, input_dict, output_dict, P_input_external_kW):
         # print("electricity_profile: ", electricity_profile)
         # P_input_external_kW = electricity_profile.iloc[:, 1].to_numpy()
         self.input_dict['P_input_external_kW'] = P_input_external_kW
-        self.power_supply_rating_MW = 15
+
+
+        # self.power_supply_rating_MW = 15 # rating of the plant powering the electrolyzer
+        self.electrolyzer_system_size_MW = 15
 
         # Uncomment line below if this model is being supplied by a 1-D np
         # array of time series external power supply (instead of reading CSV):
@@ -128,9 +131,15 @@ def system_design(self):
         system - which may consist of multiple stacks connected together in
         series, parallel, or a combination of both.
         """
-        h2_production_multiplier = (self.power_supply_rating_MW * 1000) / \
+        print("self.electrolyzer_system_size_MW: ", self.electrolyzer_system_size_MW)
+        print("self.stack_rating_kW: ", self.stack_rating_kW)
+        h2_production_multiplier = (self.electrolyzer_system_size_MW * 1000) / \
                                    self.stack_rating_kW
-        self.output_dict['electrolyzer_system_size_MW'] = math.floor(self.power_supply_rating_MW)
+        # h2_production_multiplier = self.stack_rating_kW/(self.power_supply_rating_MW * 1000)
+        # h2_production_multiplier = 1.0    
+        # print("h2_production_multiplier: ", h2_production_multiplier)      
+        # self.output_dict['electrolyzer_system_size_MW'] = math.floor(self.power_supply_rating_MW)
+        self.output_dict['electrolyzer_system_size_MW'] = self.electrolyzer_system_size_MW
         return h2_production_multiplier
 
     def cell_design(self):
@@ -381,17 +390,16 @@ def h2_production_rate(self):
         """
         # Single stack calculations:
         n_Tot = self.total_efficiency()
-
         h2_production_rate = n_Tot * ((self.N_cells *
                                        self.output_dict['current_input_external_Amps']) /
                                       (2 * self.F))  # mol/s
-
         h2_production_rate_g_s = h2_production_rate / self.moles_per_g_h2
         h2_produced_kg_hr = h2_production_rate_g_s * 3.6
 
         self.output_dict['stack_h2_produced_kg_hr'] = h2_produced_kg_hr
 
         # Total electrolyzer system calculations:
+        print("self.system_design(): ", self.system_design())
         h2_produced_kg_hr_system = self.system_design() * h2_produced_kg_hr
         self.output_dict['h2_produced_kg_hr_system'] = h2_produced_kg_hr_system
 

hybrid/hybrid_simulation.py

@@ -166,7 +166,7 @@ def calculate_installed_cost(self):
             self.solar.system_capacity_kw / 1000, storage_mw=self.storage_kw / 1000, storage_mwh=self.storage_kwh / 1000,
             storage_hours=self.storage_hours)
 
-        print("Total Solar Cost: {}, Total Wind Cost: {}, Total Storage Cost: {}".format(solar_cost, wind_cost, storage_cost))
+        # print("Total Solar Cost: {}, Total Wind Cost: {}, Total Storage Cost: {}".format(solar_cost, wind_cost, storage_cost))
         if self.solar:
             self.solar.set_total_installed_cost_dollars(solar_cost)
         if self.wind:

tools/analysis/bos/bos_lookup.py

@@ -49,7 +49,7 @@ def _lookup_costs(self, wind_mw, solar_mw, interconnection_mw):
             return 0, 0, 0
 
         search_inputs = np.array([interconnection_mw, wind_mw, solar_mw])
-        print("Search Inputs: {}".format(search_inputs))
+        # print("Search Inputs: {}".format(search_inputs))
         distance_norm = np.linalg.norm(self.contents - search_inputs, axis=1)
         min_index = np.argmin(distance_norm)
         min_distance = distance_norm[min_index]
@@ -69,9 +69,9 @@ def _lookup_costs(self, wind_mw, solar_mw, interconnection_mw):
             solar_bos_cost = vals[self.desired_output_parameters.index("Solar BOS Cost")]
 
         total_bos_cost = wind_bos_cost + solar_bos_cost
-        print("Wind BOS Cost {}, Solar BOS Cost {}, Total BOS Cost {}".format(wind_bos_cost, solar_bos_cost,
-                                                                              total_bos_cost))
-        print("MIN INDEX {}".format(min_index))
+        # print("Wind BOS Cost {}, Solar BOS Cost {}, Total BOS Cost {}".format(wind_bos_cost, solar_bos_cost,
+                                                                            #   total_bos_cost))
+        # print("MIN INDEX {}".format(min_index))
         logger.info("Total BOS Cost: {} Wind BOS Cost: {} Solar BOS Cost {}".
                     format(total_bos_cost, wind_bos_cost, solar_bos_cost))
 

tools/analysis/bos/cost_calculator.py

@@ -126,7 +126,7 @@ def calculate_total_costs(self, wind_mw, solar_mw, storage_mw=0, storage_mwh=0,
             total_solar_cost = solar_bos_cost
             total_storage_cost = storage_bos_cost
             total_project_cost = total_bos_cost
-        print("Modify costs is: {}".format(self.modify_costs))
+        # print("Modify costs is: {}".format(self.modify_costs))
         if self.modify_costs:
             logger.info('Modifying costs using selected multipliers')
             logger.info("Total Project Cost Before Modifiers: {}".format(total_project_cost))
@@ -154,7 +154,7 @@ def calculate_total_costs(self, wind_mw, solar_mw, storage_mw=0, storage_mwh=0,
             # Not modifying wind or solar costs
 
         logger.info("Total Project Cost (Installed Cost + BOS Cost): {}".format(total_project_cost))
-        print("Total Project Cost: {}".format(total_project_cost))
+        # print("Total Project Cost: {}".format(total_project_cost))
         return total_solar_cost, total_wind_cost, total_storage_cost, total_project_cost