Merge branch 'main' into updates_idaes

dispatches/case_studies/renewables_case/ConceptualDesignOptimization.ipynb

@@ -401,7 +401,8 @@
     "print(\"Cost Parameters:\")\n",
     "print(f\"Wind Capital Cost\\t\\t\\t{wind_cap_cost} $/kW\")\n",
     "print(f\"Wind Fixed Operating Cost\\t\\t{wind_op_cost} $/kW-yr\")\n",
-    "print(f\"Battery Capital Cost\\t\\t\\t{batt_cap_cost} $/kW\")\n",
+    "print(f\"Battery Power Capital Cost \\t\\t\\t{batt_cap_cost_kw} $/kW\")\n",
+    "print(f\"Battery Energy Capital Cost \\t\\t\\t{batt_cap_cost_kwh} $/kWh\")\n",
     "print(f\"PEM Capital Cost\\t\\t\\t{pem_cap_cost} $/kW\")\n",
     "print(f\"PEM Fixed Operating Cost\\t\\t{pem_op_cost} $/kW-yr\")\n",
     "print(f\"PEM Variable Operating Cost\\t\\t{pem_var_cost} $/kWh\")\n",
@@ -829,8 +830,7 @@
     "tank_type = 'simple'\n",
     "input_params['wind_resource'] = {t:\n",
     "                                    {'wind_resource_config': {\n",
-    "                                         'resource_speed': [wind_speeds[t]]\n",
-    "                                    }\n",
+    "                                        'capacity_factor': [wind_cfs[t]]}\n",
     "                                } for t in range(7)}\n",
     "                                \n",
     "mp_model = MultiPeriodModel(n_time_points=7,\n",

dispatches/case_studies/renewables_case/PEM_parametrized_bidder.py

@@ -91,7 +91,7 @@ def compute_day_ahead_bids(self, date, hour=0):
         return full_bids
 
     def compute_real_time_bids(
-        self, date, hour, _, __
+        self, date, hour, realized_day_ahead_prices, realized_day_ahead_dispatches
     ):
         """
         RT Bid: from 0 MW to (Wind Resource - PEM capacity) MW, bid $0/MWh.

dispatches/case_studies/renewables_case/README.md

@@ -0,0 +1,26 @@
+# Renewable Energy Case Study: Wind + PEM and Wind + Battery Plants
+
+This directory contains the files required for the renewable energy case studies, which use flowsheets in which any of the following 
+technologies (modeled as unit models) may be present: Wind, PV, Battery, PEM Electrolysis, Hydrogen Storage Tank and Hydrogen Turbine.
+The RE case studies focused on a Wind + PEM plant and a Wind + Battery plant, while the the industry-partnership case study looked at a PV + Battery + PEM + Hydrogen.
+Each of these three case studies has different approaches to modeling prices or the grid.
+
+## Wind + PEM case:
+
+1. Price-taker Design Optimization: `dispatches/case_studies/renewables_case/run_pricetaker_wind_PEM.py`
+2. Market Surrogate Design Optimization: `dispatches/case_studies/renewables_case/RE_surrogate_optimization_steadystate.py` 
+3. Double Loop Simulation for Validation: `dispatches/case_studies/renewables_case/run_double_loop_PEM.py`
+4. Market Surrogate Design and Validation Plotting: `dispatches/case_studies/renewables_case/SurrogateDesignResults.ipynb`
+
+## Wind + Battery case:
+1. Price-taker Design Optimization: `dispatches/case_studies/renewables_case/run_pricetaker_battery_ratio_size.py`
+2. Double Loop Simulation: `dispatches/case_studies/renewables_case/run_double_loop_battery.py`
+3. Parametrized Bidder Double Loop Simulation: `dispatches/case_studies/renewables_case/run_double_loop_battery_parametrized.py`
+
+## PV + Battery + Hydrogen case:
+1. Price-taker Design Optimization: `dispatches/case_studies/renewables_case/solar_battery_hydrogen.py`
+
+## Software and Hardware
+All the models are run on a Red Hat Enterprise Linux Server version 7.9 (Maipo). The versions for the solvers used are given below:
+- Xpress: Version 8.13.0
+- IPOPT: Version 3.13.2

dispatches/case_studies/renewables_case/RE_flowsheet.py

@@ -126,8 +126,9 @@ def add_pem(m, outlet_pressure_bar):
 
     m.fs.pem = PEM_Electrolyzer(property_package=m.fs.h2ideal_props)
 
-    # Conversion of kW to mol/sec of H2. (elec*elec_to_mol) based on H-tec design of 54.517kW-hr/kg
-    m.fs.pem.electricity_to_mol.fix(0.002527406)
+    # Conversion of kW to mol/sec of H2. (elec*elec_to_mol)
+    # 50 kWh/kg since that is the value for most NEL PEM electrolyzers: https://nelhydrogen.com/product/m-series-3/ (see M3000 data)
+    m.fs.pem.electricity_to_mol.fix(0.00275984)
     m.fs.pem.outlet.pressure.setub(max_pressure_bar * 1e5)
     m.fs.pem.outlet.pressure.fix(outlet_pressure_bar * 1e5)
     m.fs.pem.outlet.temperature.fix(pem_temp)
@@ -151,7 +152,7 @@ def add_battery(m, batt_mw):
     m.fs.battery.discharging_eta.set_value(0.95)
     m.fs.battery.dt.set_value(timestep_hrs)
     m.fs.battery.nameplate_power.fix(batt_mw * 1e3)
-    m.fs.battery.duration = Param(default=4, mutable=True, units=pyunits.kWh/pyunits.kW)
+    m.fs.battery.duration = Param(default=duration, mutable=True, units=pyunits.kWh/pyunits.kW)
     m.fs.battery.four_hr_battery = Constraint(expr=m.fs.battery.nameplate_power * m.fs.battery.duration == m.fs.battery.nameplate_energy)
     return m.fs.battery
 

dispatches/case_studies/renewables_case/RE_surrogate_optimization_steadystate.py

@@ -45,6 +45,11 @@
                                 initialize_mp, wind_battery_pem_model, wind_battery_pem_mp_block
 
 
+# RT market only or Both RT and DA markets
+rt_market_only = True
+include_wind_capital_cost = False
+shortfall = 1000
+
 # path for folder that has surrogate models
 re_nn_dir = Path(__file__).parent / "data" / "steady_state_surrogate"
 
@@ -57,16 +62,21 @@ def load_surrogate_model(re_nn_dir):
     pem_clusters_mean = centers[:, 1]
     resource_clusters_mean = centers[:, 2]
 
-    with open(re_nn_dir / "revenue" / "RE_revenue_params_2_25.json", 'r') as f:
-        rev_data = json.load(f)
-
     # load keras neural networks
     # Input variables are PEM bid price, PEM MW, Reserve Factor and Load Shed Price
-    nn_rev = keras.models.load_model(re_nn_dir / "revenue" / "RE_revenue_2_25")
-    nn_dispatch = keras.models.load_model(re_nn_dir / "dispatch_frequency" / "ss_surrogate_model_wind_pmax")
-
     with open(re_nn_dir / "dispatch_frequency" / "ss_surrogate_param_wind_pmax.json", 'r') as f:
         dispatch_data = json.load(f)
+    nn_dispatch = keras.models.load_model(re_nn_dir / "dispatch_frequency" / "ss_surrogate_model_wind_pmax")
+
+    if rt_market_only:
+        rev_data_f = re_nn_dir / "rt_revenue" / "RE_RT_revenue_params_2_25.json"
+        nn_rev = keras.models.load_model(re_nn_dir / "rt_revenue" / "RE_RT_revenue_2_25")
+    else:
+        rev_data_f = re_nn_dir / "revenue" / "RE_revenue_params_2_25.json"
+        nn_rev = keras.models.load_model(re_nn_dir / "revenue" / "RE_revenue_2_25")
+
+    with open(rev_data_f, 'rb') as f:
+        rev_data = json.load(f)
 
     # load keras models and create OMLT NetworkDefinition objects
     #Revenue model definition
@@ -97,10 +107,10 @@ def conceptual_design_dynamic_RE(input_params, PEM_bid=None, PEM_MW=None, verbos
     # add surrogate input to the model
     m.wind_system_capacity = Var(domain=NonNegativeReals, bounds=(100 * 1e3, 1000 * 1e3), initialize=input_params['wind_mw'] * 1e3)
 
-    m.pem_system_capacity = Var(domain=NonNegativeReals, bounds=(127.05 * 1e3, 423.5 * 1e3), initialize=input_params['pem_mw'] * 1e3, units=pyunits.kW)
+    m.pem_system_capacity = Var(domain=NonNegativeReals, bounds=(127.5 * 1e3, 423.5 * 1e3), initialize=input_params['pem_mw'] * 1e3, units=pyunits.kW)
     m.pem_bid = Var(within=NonNegativeReals, bounds=(15, 45), initialize=45)                    # Energy Bid $/MWh
     m.reserve_percent = Param(within=NonNegativeReals, initialize=15)   # Reserves Fraction on Grid
-    m.shortfall_price = Param(within=NonNegativeReals, initialize=1000)     # Energy price during load shed
+    m.shortfall_price = Param(within=NonNegativeReals, initialize=shortfall)     # Energy price during load shed
 
     inputs = [m.pem_bid, m.pem_system_capacity * 1e-3 / 847 * m.wind_system_capacity * 1e-3, m.reserve_percent, m.shortfall_price]
 
@@ -213,7 +223,7 @@ def conceptual_design_dynamic_RE(input_params, PEM_bid=None, PEM_MW=None, verbos
         scenario_models.append(scenario_model)
 
     m.plant_cap_cost = Expression(
-        expr=input_params['wind_cap_cost'] * m.wind_system_capacity + input_params['pem_cap_cost'] * m.pem_system_capacity)
+        expr=input_params['wind_cap_cost'] * m.wind_system_capacity * int(include_wind_capital_cost) + input_params['pem_cap_cost'] * m.pem_system_capacity)
     m.annual_fixed_cost = pyo.Expression(
         expr=m.wind_system_capacity * input_params["wind_op_cost"] + m.pem_system_capacity * input_params["pem_op_cost"])
     m.plant_operation_cost = Expression(
@@ -222,7 +232,8 @@ def conceptual_design_dynamic_RE(input_params, PEM_bid=None, PEM_MW=None, verbos
         expr=sum(scenario_models[i].hydrogen_revenue for i in range(num_rep_days)))
 
     m.NPV = Expression(expr=-m.plant_cap_cost + PA * (m.rev + m.hydrogen_rev - m.plant_operation_cost - m.annual_fixed_cost))
-    m.obj = Objective(expr=-m.NPV * 1e-8)
+    m.NPV_ann = Expression(expr=-m.plant_cap_cost / PA + (m.rev + m.hydrogen_rev - m.plant_operation_cost - m.annual_fixed_cost))
+    m.obj = Objective(expr=-m.NPV * 1e-7)
 
     return m, num_rep_days
 
@@ -250,7 +261,8 @@ def record_result(m, num_rep_days):
         "pem_bid": value(m.pem_bid),
         "e_revenue": value(m.rev),
         "h_revenue": value(m.hydrogen_rev),
-        "NPV": value(m.NPV)
+        "NPV": value(m.NPV),
+        "NPV_ann": value(m.NPV_ann)
     }
 
     for day in range(num_rep_days):
@@ -263,6 +275,7 @@ def record_result(m, num_rep_days):
     print("Plant Hydrogen Revenue Annual = ${}".format(value(m.hydrogen_rev)))
     print("Plant Total Revenue Annual = ${}".format(value(m.rev + m.hydrogen_rev)))
     print("Plant NPV = ${}".format(value(m.NPV)))
+    print("Plant NPV Annualized = ${}".format(value(m.NPV_ann)))
 
     print('----------')
     for i in range(num_rep_days):
@@ -279,8 +292,6 @@ def run_design(PEM_bid=None, PEM_size=None):
     model, n_rep_days = conceptual_design_dynamic_RE(default_input_params, PEM_bid=PEM_bid, PEM_MW=PEM_size, verbose=False)
     nlp_solver = SolverFactory('ipopt')
     nlp_solver.options['max_iter'] = 8000
-    nlp_solver.options['acceptable_tol'] = 1e-8
-    nlp_solver.options['bound_push'] = 1e-9
     res = nlp_solver.solve(model, tee=True)
     if res.Solver.status != 'ok':
         solve_log = idaeslog.getInitLogger("infeasibility", idaeslog.INFO, tag="properties")
@@ -314,7 +325,7 @@ def run_design(PEM_bid=None, PEM_size=None):
     "design_opt": True,
     "extant_wind": True,        # fixed because parameter sweeps didn't change wind size
 
-    "wind_cap_cost": wind_cap_cost,
+    "wind_cap_cost": wind_cap_cost if include_wind_capital_cost else 0,
     "wind_op_cost": wind_op_cost,
     "pem_cap_cost": pem_cap_cost,
     "pem_op_cost": pem_op_cost,
@@ -323,8 +334,8 @@ def run_design(PEM_bid=None, PEM_size=None):
 
 
 if __name__ == "__main__":
-    # result = run_design()
-    # exit()
+    result = run_design()
+    exit()
 
     import multiprocessing as mp
     from itertools import product
@@ -333,8 +344,8 @@ def run_design(PEM_bid=None, PEM_size=None):
     sizes = np.linspace(127.05, 423.5, 15)
     inputs = product(bids, sizes)
 
-    with mp.Pool(processes=4) as p:
+    with mp.Pool(processes=24) as p:
         res = p.starmap(run_design, inputs)
 
     df = pd.DataFrame(res)
-    df.to_csv("surrogate_results_ss.csv")
+    df.to_csv(f"wind_PEM/surrogate_results_ss_rt_{shortfall}.csv")

dispatches/case_studies/renewables_case/battery_parametrized_bidder.py

@@ -0,0 +1,124 @@
+#################################################################################
+# DISPATCHES was produced under the DOE Design Integration and Synthesis Platform
+# to Advance Tightly Coupled Hybrid Energy Systems program (DISPATCHES), and is
+# copyright (c) 2020-2023 by the software owners: The Regents of the University
+# of California, through Lawrence Berkeley National Laboratory, National
+# Technology & Engineering Solutions of Sandia, LLC, Alliance for Sustainable
+# Energy, LLC, Battelle Energy Alliance, LLC, University of Notre Dame du Lac, et
+# al. All rights reserved.
+#
+# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and license
+# information, respectively. Both files are also available online at the URL:
+# "https://github.com/gmlc-dispatches/dispatches".
+#################################################################################
+import numpy as np
+from idaes.apps.grid_integration.bidder import convert_marginal_costs_to_actual_costs, tx_utils
+from dispatches.workflow.parametrized_bidder import ParametrizedBidder
+
+
+class FixedParametrizedBidder(ParametrizedBidder):
+
+    """
+    Template class for bidders that use fixed parameters. 
+    The functions for computing the day ahead and real time bids do not use any information from Prescient,
+    and only depend on internal system information such as wind capacity factors, and on bid parameters.
+    """
+
+    def __init__(
+        self,
+        bidding_model_object,
+        day_ahead_horizon,
+        real_time_horizon,
+        solver,
+        forecaster,
+        storage_marginal_cost,
+        storage_mw
+    ):
+        super().__init__(bidding_model_object,
+                         day_ahead_horizon,
+                         real_time_horizon,
+                         solver,
+                         forecaster)
+        self.wind_marginal_cost = 0
+        self.wind_mw = self.bidding_model_object._wind_pmax_mw
+        self.storage_marginal_cost = storage_marginal_cost
+        self.storage_mw = storage_mw
+
+    def compute_day_ahead_bids(self, date, hour=0):
+        """
+        Day ahead bid has two parts: 
+            1. the part of the DA wind energy that is able to be stored in the battery is at the higher "storage_marginal_cost"
+            2. the remainder of the wind energy is bid at the wind marginal cost of $0/MWh
+
+        For each time period in the day ahead horizon, the marginal cost bid is assembled as the two parts. 
+        Then the marginal costs are converted to actual costs, as expected by Prescient.
+        """
+        gen = self.generator
+        forecast = self.forecaster.forecast_day_ahead_capacity_factor(date, hour, gen, self.day_ahead_horizon)
+
+        full_bids = {}
+
+        for t_idx in range(self.day_ahead_horizon):
+            da_wind = forecast[t_idx] * self.wind_mw
+            p_max = max(da_wind, self.storage_mw)
+            bids = [(0, 0), (max(0, da_wind - self.storage_mw), 0), (p_max, self.storage_marginal_cost)]
+            cost_curve = convert_marginal_costs_to_actual_costs(bids)
+
+            temp_curve = {
+                    "data_type": "cost_curve",
+                    "cost_curve_type": "piecewise",
+                    "values": cost_curve,
+            }
+            tx_utils.validate_and_clean_cost_curve(
+                curve=temp_curve,
+                curve_type="cost_curve",
+                p_min=0,
+                p_max=max([p[0] for p in cost_curve]),
+                gen_name=gen,
+                t=t_idx,
+            )
+
+            t = t_idx + hour
+            full_bids[t] = {}
+            full_bids[t][gen] = {}
+            full_bids[t][gen]["p_cost"] = cost_curve
+            full_bids[t][gen]["p_min"] = 0
+            full_bids[t][gen]["p_max"] = p_max
+            full_bids[t][gen]["startup_capacity"] = p_max
+            full_bids[t][gen]["shutdown_capacity"] = p_max
+
+        self._record_bids(full_bids, date, hour, Market="Day-ahead")
+        return full_bids
+
+    def compute_real_time_bids(
+        self, date, hour, realized_day_ahead_prices, realized_day_ahead_dispatches
+    ):
+        """
+        Real time bid has two parts: 
+            1. the part of the RT wind energy that is able to be stored in the battery is at the higher "storage_marginal_cost"
+            2. the remainder of the wind energy is bid at the wind marginal cost of $0/MWh
+
+        For each time period in the day ahead horizon, the marginal cost bid is assembled as the two parts. 
+        Then the marginal costs are converted to actual costs, as expected by Prescient.
+        """
+        gen = self.generator
+        forecast = self.forecaster.forecast_real_time_capacity_factor(date, hour, gen, self.day_ahead_horizon)
+
+        full_bids = {}
+
+        for t_idx in range(self.real_time_horizon):
+            rt_wind = forecast[t_idx] * self.wind_mw
+            p_max = max(rt_wind, self.storage_mw)
+            bids = [(0, 0),  (max(0, rt_wind - self.storage_mw), 0), (p_max, self.storage_marginal_cost)]
+
+            t = t_idx + hour
+            full_bids[t] = {}
+            full_bids[t][gen] = {}
+            full_bids[t][gen]["p_cost"] = convert_marginal_costs_to_actual_costs(bids)
+            full_bids[t][gen]["p_min"] = 0
+            full_bids[t][gen]["p_max"] = p_max
+            full_bids[t][gen]["startup_capacity"] = p_max
+            full_bids[t][gen]["shutdown_capacity"] = p_max
+
+        self._record_bids(full_bids, date, hour, Market="Real-time")
+        return full_bids