Merge pull request #288 from jaredthomas68/greenheart

Greenheart

greenheart/simulation/technologies/hydrogen/electrolysis/PEM_costs_Singlitico_model.py

@@ -66,6 +66,7 @@ def run(
             tuple: CapEx and OpEx costs for a single electrolyzer.
         """
         capex = self.calc_capex(P_elec, RC_elec)
+
         opex = self.calc_opex(P_elec, capex)
 
         return capex, opex

greenheart/tools/eco/hopp_mgmt.py

@@ -18,7 +18,11 @@ def setup_hopp(
     show_plots=False,
     save_plots=False,
 ):
-    hopp_site = SiteInfo(**hopp_config["site"], desired_schedule=[greenheart_config["electrolyzer"]["rating"]]*8760)
+    if "battery" in hopp_config["technologies"].keys() and \
+        ("desired_schedule" not in hopp_config["site"].keys() or hopp_config["site"]["desired_schedule"] == []):
+        hopp_config["site"]["desired_schedule"] = [greenheart_config["electrolyzer"]["rating"]]*8760
+    # hopp_site = SiteInfo(**hopp_config["site"], desired_schedule=[greenheart_config["electrolyzer"]["rating"]]*8760)
+    hopp_site = SiteInfo(**hopp_config["site"])
 
     # adjust mean wind speed if desired
     wind_data = hopp_site.wind_resource._data['data']
@@ -123,7 +127,7 @@ def run_hopp(hopp_config, hopp_site, project_lifetime, verbose=False):
 
     if "wave" in hi.system.technologies.keys():
         hi.system.wave.create_mhk_cost_calculator(wave_cost_dict)
-
+        
     hi.simulate(project_life=project_lifetime)
 
     # store results for later use

greenheart/tools/optimization/openmdao.py

@@ -5,27 +5,54 @@
 from electrolyzer import run_lcoh
 
 from shapely.geometry import Polygon, Point
+from hopp.simulation import HoppInterface
 
 class HOPPComponent(om.ExplicitComponent):
 
     def initialize(self):
-        self.options.declare("hi", recordable=False)
-        self.options.declare("turbine_x_init")
-        self.options.declare("turbine_y_init")
-        self.options.declare("verbose")
+        self.options.declare("hi", types=HoppInterface, recordable=False, desc="HOPP Interface class instance")
+        self.options.declare("turbine_x_init", types=(list, type(np.array(0))), desc="Initial turbine easting locations in m")
+        self.options.declare("turbine_y_init", types=(list, type(np.array(0))), desc="Initial turbine northing locations in m")
+        self.options.declare("verbose", default=False, types=bool, desc="Whether or not to print a bunch of stuff")
+        self.options.declare("design_variables", 
+                            #  values=["turbine_x", "turbine_y", "pv_capacity_kw", "wind_rating_kw", "electrolyzer_rating_kw", "battery_capacity_kw", "battery_capacity_kwh"], 
+                             types=list, 
+                             desc="List of design variables that should be included",
+                             default=["turbine_x", "turbine_y"],
+                             recordable=False)
 
     def setup(self):
 
-        self.add_input("turbine_x", val=self.options["turbine_x_init"], units="m")
-        self.add_input("turbine_y", val=self.options["turbine_y_init"], units="m")
-        # self.add_input("battery_capacity_kw", val=15000, units="kW")
-        # self.add_input("battery_capacity_kwh", val=15000, units="kW*h")
-        # self.add_input("electrolyzer_rating_kw", val=15000, units="kW*h")
-        # self.add_input("pv_rating_kw", val=15000, units="kW")
-        # self.add_input("wind_rating_kw", val=150000, units="kW")
+        if "turbine_x" in self.options["design_variables"]:
+           self.add_input("turbine_x", val=self.options["turbine_x_init"], units="m")
+        if "turbine_y" in self.options["design_variables"]:
+           self.add_input("turbine_y", val=self.options["turbine_y_init"], units="m")
+        if "wind_rating_kw" in self.options["design_variables"]:
+            self.add_input("wind_rating_kw", val=150000, units="kW")
+        if "pv_capacity_kw" in self.options["design_variables"]:
+            self.add_input("pv_capacity_kw", val=15000, units="kW")
+        if "battery_capacity_kw" in self.options["design_variables"]:
+            self.add_input("battery_capacity_kw", val=15000, units="kW")
+        if "battery_capacity_kwh" in self.options["design_variables"]:
+            self.add_input("battery_capacity_kwh", val=15000, units="kW*h")
+
+        technologies = self.options["hi"].configuration["technologies"]
+
+        if "pv" in technologies.keys():
+            self.add_output("pv_capex", val=0.0, units="USD")
+            # self.add_output("pv_opex", val=0.0, units="USD")
+        if "wind" in technologies.keys():
+            self.add_output("wind_capex", val=0.0, units="USD")
+            # self.add_output("wind_opex", val=0.0, units="USD")
+        if "battery" in technologies.keys():
+            self.add_output("battery_capex", val=0.0, units="USD")
+            # self.add_output("battery_opex", val=0.0, units="USD")
+
+        self.add_output("hybrid_electrical_generation_capex", units="USD")
+        self.add_output("hybrid_electrical_generation_opex", units="USD")
         self.add_output("aep", units="kW*h")
         self.add_output("lcoe_real", units="USD/(MW*h)")
-        self.add_output("p_wind", units="kW", val=np.zeros(8760))
+        self.add_output("power_signal", units="kW", val=np.zeros(8760))
 
     def compute(self, inputs, outputs):
 
@@ -36,21 +63,55 @@ def compute(self, inputs, outputs):
 
         hi = self.options["hi"]
 
-        hi.system.wind._system_model.fi.reinitialize(layout_x=inputs["turbine_x"], layout_y=inputs["turbine_y"], time_series=True)
-
+        if any(x in ["wind_rating_kw", "pv_capacity_kw", "battery_capacity_kw", "battery_capacity_kwh"] for x in inputs):
+            technologies = hi.configuration["technologies"]
+            if "wind_rating_kw" in inputs:
+                raise(NotImplementedError("wind_rating_kw has not be fully implemented as a design variable"))
+            if "pv_capacity_kw" in inputs:
+                technologies["pv"]["system_capacity_kw"] = float(inputs["pv_capacity_kw"])
+            if "battery_capacity_kw" in inputs:
+                technologies["battery"]["system_capacity_kw"] = float(inputs["battery_capacity_kw"])
+            if "battery_capacity_kwh" in inputs:
+                technologies["battery"]["system_capacity_kwh"] = float(inputs["battery_capacity_kwh"])
+
+            configuration = hi.configuration
+            configuration["technologies"] = technologies
+            hi.reinitialize(configuration)
+
+        if ("turbine_x" in inputs) or ("turbine_y" in inputs):
+            if "turbine_x" not in inputs:
+                hi.system.wind._system_model.fi.reinitialize(layout_y=inputs["turbine_y"], time_series=True)
+            elif "turbine_y" not in inputs:
+                hi.system.wind._system_model.fi.reinitialize(layout_x=inputs["turbine_x"], time_series=True)
+            else:
+                hi.system.wind._system_model.fi.reinitialize(layout_x=inputs["turbine_x"], layout_y=inputs["turbine_y"], time_series=True)
+
         # run simulation
         hi.simulate(25)
 
         if self.options["verbose"]:
             print(f"obj: {hi.system.annual_energies.hybrid}")
 
-        # get result
+        # get results
+        if "pv" in technologies.keys():
+            outputs["pv_capex"] = hi.system.pv.total_installed_cost
+            # outputs["pv_opex"] = hi.system.pv.om_total_expense[1]
+        if "wind" in technologies.keys():
+            outputs["wind_capex"] = hi.system.wind.total_installed_cost
+            # outputs["wind_opex"] = hi.system.wind.om_total_expense[1]
+        if "battery" in technologies.keys():
+            outputs["battery_capex"] = hi.system.battery.total_installed_cost
+            # outputs["battery_opex"] = hi.system.battery.om_total_expense[1]
+
+        outputs["hybrid_electrical_generation_capex"] = hi.system.cost_installed["hybrid"]
+        # outputs["hybrid_electrical_generation_opex"] = hi.system.om_total_expenses[1]
+
         outputs["aep"] = hi.system.annual_energies.hybrid
         outputs["lcoe_real"] = hi.system.lcoe_real.hybrid/100.0 # convert from cents/kWh to USD/kWh
-        outputs["p_wind"] = hi.system.grid.generation_profile[0:8760]
+        outputs["power_signal"] = hi.system.grid.generation_profile[0:8760]
 
     def setup_partials(self):
-        self.declare_partials('*', '*', method='fd', form="forward")
+        self.declare_partials(['lcoe_real', 'power_signal'], '*', method='fd', form="forward")
 
 class TurbineDistanceComponent(om.ExplicitComponent):
 
@@ -99,7 +160,7 @@ def compute(self, inputs, outputs):
         hi = self.options["hopp_interface"]
 
         # get polygon for boundary
-        boundary_polygon = Polygon(hi.system.site.vertices)     
+        boundary_polygon = Polygon(hi.system.site.vertices) 
         # check if turbines are inside polygon and get distance
         for i in range(0, self.n_distances):
             point = Point(inputs["turbine_x"][i], inputs["turbine_y"][i])
@@ -122,25 +183,38 @@ def initialize(self):
         self.options.declare("h2_modeling_options")
         self.options.declare("h2_opt_options")
         self.options.declare("modeling_options")
+        self.options.declare("design_variables", 
+                             # ["electrolyzer_rating_kw"], 
+                             types=list, 
+                             desc="List of design variables that should be included",
+                             default=[],
+                             recordable=False)
 
     def setup(self):
-        self.add_input("p_wind", val=np.zeros(8760), units="W")
+        self.add_input("power_signal", val=np.zeros(8760), units="W")
         self.add_input("lcoe_real", units="USD/kW/h")
+
+        if "electrolyzer_rating_kw" in self.options["design_variables"]:
+            self.add_input("electrolyzer_rating_kw", val=15000, units="kW*h")
+
         if self.options["h2_opt_options"]["control"]["system_rating_MW"]["flag"] \
             or self.options["modeling_options"]["rating_equals_turbine_rating"]:
             self.add_input("system_rating_MW", units="MW", val=self.options["h2_modeling_options"]["electrolyzer"]["control"]["system_rating_MW"])
         self.add_output("h2_produced", units="kg")
         self.add_output("max_curr_density", units="A/cm**2")
-        self.add_output("capex", units="USD")
-        self.add_output("opex", units="USD")
-        self.add_output("lcoh", units="USD/kg")
+        self.add_output("electrolyzer_capex", units="USD")
+        self.add_output("electrolyzer_opex", units="USD")
+        self.add_output("lcoh", units="USD/kg") 
 
     def compute(self, inputs, outputs):
         # Set electrolyzer parameters from model inputs
-        power_signal = inputs["p_wind"]
+        power_signal = inputs["power_signal"]
         lcoe_real = inputs["lcoe_real"][0]
 
-        if self.options["h2_opt_options"]["control"]["system_rating_MW"]["flag"] \
+        if "electrolyzer_rating_kw" in inputs:
+            self.options["h2_modeling_options"]["electrolyzer"]["control"]["system_rating_MW"] = inputs["electrolyzer_rating_kw"]
+
+        elif self.options["h2_opt_options"]["control"]["system_rating_MW"]["flag"] \
             or self.options["modeling_options"]["rating_equals_turbine_rating"]:
             self.options["h2_modeling_options"]["electrolyzer"]["control"]["system_rating_MW"] = inputs["system_rating_MW"][0]
 
@@ -170,9 +244,9 @@ def compute(self, inputs, outputs):
 
         outputs["h2_produced"] = h2_prod
         outputs["max_curr_density"] = max_curr_density
-        outputs["capex"] = capex
-        outputs["opex"] = opex
+        outputs["electrolyzer_capex"] = capex
+        outputs["electrolyzer_opex"] = opex
         outputs["lcoh"] = lcoh
 
     def setup_partials(self):
-        self.declare_partials('*', '*', method='fd', form='forward')
+        self.declare_partials('lcoh', '*', method='fd', form='forward')

hopp/simulation/hopp_interface.py

@@ -2,7 +2,7 @@
 from pathlib import Path
 from typing import Union, TYPE_CHECKING
 
-from hopp.simulation.hopp import Hopp
+from hopp.simulation.hopp import Hopp, SiteInfo
 
 # avoid potential circular dep
 if TYPE_CHECKING:
@@ -42,8 +42,8 @@ def __init__(self, configuration: Union[dict, str, Path]):
         elif isinstance(self.configuration, dict):
             self.hopp = Hopp.from_dict(self.configuration)
 
-    def reinitialize(self):
-        pass
+    def reinitialize(self, configuration: Union[dict, str, Path]):
+        self.__init__(configuration)
 
     def simulate(self, project_life: int = 25, lifetime_sim: bool = False):
         self.hopp.simulate(project_life, lifetime_sim)

tests/greenheart/inputs/hopp_config.yaml

@@ -11,7 +11,6 @@ site: #!include flatirons_site.yaml
         - [600.0, 600.0]
         - [600.0, 0.0]
     urdb_label: "5ca4d1175457a39b23b3d45e"
-
   solar_resource_file: ""
   wind_resource_file: "" #
   wave_resource_file: ""
@@ -25,9 +24,9 @@ site: #!include flatirons_site.yaml
   wind_resource_origin: "WTK"
 
 technologies:
-  # pv:
-  #   system_capacity_kw: 50000 #
-  #   # dc_degradation: [0] * 25 #
+  pv:
+    system_capacity_kw: 50000 #
+    # dc_degradation: [0] * 25 #
   wind:
     num_turbines: 6 #
     turbine_rating_kw: 5000.0 #

tests/greenheart/test_hydrogen/input_files/floris/floris_input_osw_18MW.yaml

@@ -25,6 +25,46 @@ flow_field:
   - 8.0
   wind_veer: 0.0
 
+###
+# Configure the turbine types and their placement within the wind farm.
+farm:
+
+  ###
+  # Coordinates for the turbine locations in the x-direction which is typically considered
+  # to be the streamwise direction (left, right) when the wind is out of the west.
+  # The order of the coordinates here corresponds to the index of the turbine in the primary
+  # data structures.
+  layout_x:
+  - 0.0
+  - 630.0
+  - 1260.0
+  - 2000.0
+  - 2200.0
+  - 2400.0
+
+  ###
+  # Coordinates for the turbine locations in the y-direction which is typically considered
+  # to be the spanwise direction (up, down) when the wind is out of the west.
+  # The order of the coordinates here corresponds to the index of the turbine in the primary
+  # data structures.
+  layout_y:
+  - 0.0
+  - 0.0
+  - 0.0
+  - 0.0
+  - 0.0
+  - 0.0
+
+  ###
+  # Listing of turbine types for placement at the x and y coordinates given above.
+  # The list length must be 1 or the same as ``layout_x`` and ``layout_y``. If it is a
+  # single value, all turbines are of the same type. Otherwise, the turbine type
+  # is mapped to the location at the same index in ``layout_x`` and ``layout_y``.
+  # The types can be either a name included in the turbine_library or
+  # a full definition of a wind turbine directly.
+  turbine_type:
+  - nrel_5MW
+
 wake:
   model_strings:
     combination_model: sosfs # TODO followup on why not using linear free-stream superposition? - ask Chris Bay and Gen S.

tests/greenheart/test_hydrogen/input_files/plant/hopp_config_wind_wave_solar_battery.yaml

@@ -88,7 +88,7 @@ config:
     wind:
       skip_financial: true
   dispatch_options:
-    battery_dispatch: heuristic #heuristic
+    battery_dispatch: load_following_heuristic
     solver: cbc
     n_look_ahead_periods: 48
     grid_charging: false

tests/greenheart/test_hydrogen/test_greenheart_system.py

@@ -130,9 +130,9 @@ def setUpClass(self):
                                     output_level=4)
 
     def test_lcoh(self):
-        assert self.lcoh == approx(13.22669818008385) #TODO base this test value on something. Currently just based on output at writing.
+        assert self.lcoh == approx(13.240698497098025) #TODO base this test value on something. Currently just based on output at writing.
     def test_lcoe(self):
-        assert self.lcoe == approx(0.13955940183722207) # TODO base this test value on something. Currently just based on output at writing.
+        assert self.lcoe == approx(0.1401530976083388) # TODO base this test value on something. Currently just based on output at writing.
 
 # class TestSimulationWindOnshore(unittest.TestCase):
 #     def setUp(self) -> None: