39 add integration test for longer timescale - checking for merge conflicts

electrolyzer/cell.py

@@ -16,9 +16,9 @@
 
 def electrolyzer_model(X, a, b, c, d, e, f):
     """
-    Given a power input, temperature, and set of coefficients, returns current.
-    Coefficients can be determined using non-linear least squares fit (see
-    `ElectrolyzerCell.create_polarization`).
+    Given a power input (kW), temperature (C), and set of coefficients, returns
+    current (A).  Coefficients can be determined using non-linear least squares
+    fit (see `Stack.create_polarization`).
     """
     P, T = X
     I = a * (P**2) + b * T**2 + c * P * T + d * P + e * T + f

electrolyzer/glue_code/optimization.py

@@ -0,0 +1,103 @@
+import copy
+
+from scipy.optimize import fsolve
+
+from electrolyzer import Stack
+
+
+def calc_rated_system(modeling_options: dict):
+    """
+    Calculates number of stacks and stack power rating (kW) to match a desired
+    system rating (MW).
+
+    Args:
+        modeling_options (dict): An options Dict compatible with the modeling schema
+    """
+    options = copy.deepcopy(modeling_options)
+
+    system_rating_kW = options["electrolyzer"]["control"]["system_rating_MW"] * 1e3
+    stack_rating_kW = options["electrolyzer"]["stack"]["stack_rating_kW"]
+
+    # determine number of stacks (int) closest to stack rating (float)
+    n_stacks = round(system_rating_kW / stack_rating_kW)
+    options["electrolyzer"]["control"]["n_stacks"] = n_stacks
+
+    # determine new desired rating to adjust parameters for
+    new_rating = system_rating_kW / n_stacks
+    options["electrolyzer"]["stack"]["stack_rating_kW"] = new_rating
+
+    # solve for new stack rating (modifies dict)
+    calc_rated_stack(options)
+
+    return options
+
+
+def _solve_rated_stack(desired_rating: float, stack: Stack):
+    cell_area = stack.cell.cell_area
+
+    # root finding function
+    def calc_rated_power_diff(cell_area: float):
+        stack.cell.cell_area = cell_area
+        p_rated = stack.calc_stack_power(stack.max_current)
+
+        return p_rated - desired_rating
+
+    return fsolve(calc_rated_power_diff, cell_area)
+
+
+def calc_rated_stack(modeling_options: dict):
+    """
+    For a given model specification, determines a configuration that meets the
+    desired stack rating (kW). Only modifies `n_cells` and `cell_area`.
+
+    NOTE: This is a naive approach: it is only concerned with achieving the desired
+        power rating. Any checks on the validity of the resulting design must be
+        performed by the user.
+
+    Args:
+        modeling_options (dict): An options Dict compatible with the modeling schema
+    """
+    options = modeling_options["electrolyzer"]["stack"]
+    options["dt"] = modeling_options["electrolyzer"]["dt"]
+    stack = Stack.from_dict(options)
+
+    n_cells = stack.n_cells
+
+    # start with an initial calculation of stack power to compare with desired
+    stack_p = stack.calc_stack_power(stack.max_current)
+    desired_rating = stack.stack_rating_kW
+
+    stack_p_prev = stack_p
+    # nudge cell count up or down until it overshoots
+    if stack_p > desired_rating:
+        while stack_p > desired_rating:
+            n_cells -= 1
+            stack.n_cells = n_cells
+            stack_p_prev = stack_p
+            stack_p = stack.calc_stack_power(stack.max_current)
+
+        # choose n_cells with closest resulting power rating
+        if abs(desired_rating - stack_p_prev) < abs(desired_rating - stack_p):
+            stack.n_cells += 1
+
+    elif stack_p < desired_rating:
+        while stack_p < desired_rating:
+            n_cells += 1
+            stack.n_cells = n_cells
+            stack_p_prev = stack_p
+            stack_p = stack.calc_stack_power(stack.max_current)
+
+        # choose n_cells with closest resulting power rating
+        if abs(desired_rating - stack_p_prev) < abs(desired_rating - stack_p):
+            stack.n_cells -= 1
+
+    # solve for optimal stack
+    cell_area = _solve_rated_stack(desired_rating, stack)
+
+    # recalc stack power
+    stack.cell.cell_area = cell_area[0]
+    stack_p = stack.calc_stack_power(stack.max_current)
+
+    modeling_options["electrolyzer"]["stack"]["cell_area"] = cell_area[0]
+    modeling_options["electrolyzer"]["stack"]["n_cells"] = n_cells
+    modeling_options["electrolyzer"]["stack"]["stack_rating_kW"] = stack_p

electrolyzer/glue_code/run_electrolyzer.py

@@ -9,37 +9,10 @@
 import electrolyzer.inputs.validation as val
 from electrolyzer import Supervisor
 
+from .optimization import calc_rated_system
 
-def run_electrolyzer(input_modeling, power_signal):
-    """
-    Runs an electrolyzer simulation based on a YAML configuration file and power
-    signal input.
-
-    Args:
-        input_modeling (`str` or `dict`): filepath specifying the YAML config
-            file, OR a dict representing a validated YAML config.
-        power_signal (`list`): An array representing power input
-
-    Returns:
-        `Supervisor`: the instance used to run the simulation
-        `pandas.DataFrame`: a `DataFrame` representing the time series output
-            for the system, including values for each electrolyzer stack
-    """
-    err_msg = "Model input must be a str or dict object"
-    assert isinstance(
-        input_modeling,
-        (
-            str,
-            dict,
-        ),
-    ), err_msg
-
-    if isinstance(input_modeling, str):
-        # Parse/validate yaml configuration
-        modeling_options = val.load_modeling_yaml(input_modeling)
-    else:
-        modeling_options = input_modeling
 
+def _run_electrolyzer_full(modeling_options, power_signal):
     # Initialize system
     elec_sys = Supervisor.from_dict(modeling_options["electrolyzer"])
 
@@ -50,10 +23,12 @@ def run_electrolyzer(input_modeling, power_signal):
     tot_kg = np.zeros((len(power_signal)))
     cycles = np.zeros((elec_sys.n_stacks, len(power_signal)))
     uptime = np.zeros((elec_sys.n_stacks, len(power_signal)))
+    current_density = np.zeros((elec_sys.n_stacks, len(power_signal)))
     p_in = []
 
     # Run electrolyzer simulation
     for i in range(len(power_signal)):
+        # TODO: replace with proper logging
         # if (i % 1000) == 0:
         #     print('Progress', i)
         # print(i)
@@ -62,27 +37,109 @@ def run_electrolyzer(input_modeling, power_signal):
         )
         p_in.append(power_signal[i] / elec_sys.n_stacks / 1000)
 
-        tot_kg[i] = np.copy(loop_H2)
-        curtailment[i] = np.copy(curtailed) / 1000000
+        tot_kg[i] = loop_H2
+        curtailment[i] = curtailed / 1000000
         for j in range(elec_sys.n_stacks):
+            stack = elec_sys.stacks[j]
             kg_rate[j, i] = loop_h2_mfr[j]
-            degradation[j, i] = elec_sys.stacks[j].V_degradation
-            cycles[j, i] = elec_sys.stacks[j].cycle_count
-            uptime[j, i] = elec_sys.stacks[j].uptime
+            degradation[j, i] = stack.V_degradation
+            cycles[j, i] = stack.cycle_count
+            uptime[j, i] = stack.uptime
+            current_density[j, i] = stack.I / stack.cell.cell_area
 
     # Collect results into a DataFrame
-    results_df = pd.DataFrame(index=range(len(power_signal)))
+    results_df = pd.DataFrame(
+        {
+            "power_signal": power_signal,
+            "curtailment": curtailment,
+            "kg_rate": tot_kg,
+        }
+    )
 
-    results_df["power_signal"] = power_signal
-    results_df["curtailment"] = curtailment
-    results_df["kg_rate"] = tot_kg
+    # for efficiency reasons, create a df for each stack, then concat all at the end
+    stack_dfs = []
 
     for i, stack in enumerate(elec_sys.stacks):
         id = i + 1
-        results_df[f"stack_{id}_deg"] = degradation[i, :]
-        results_df[f"stack_{id}_fatigue"] = stack.fatigue_history
-        results_df[f"stack_{id}_cycles"] = cycles[i, :]
-        results_df[f"stack_{id}_uptime"] = uptime[i, :]
-        results_df[f"stack_{id}_kg_rate"] = kg_rate[i, :]
+        stack_df = pd.DataFrame(
+            {
+                f"stack_{id}_deg": degradation[i, :],
+                f"stack_{id}_fatigue": stack.fatigue_history,
+                f"stack_{id}_cycles": cycles[i, :],
+                f"stack_{id}_uptime": uptime[i, :],
+                f"stack_{id}_kg_rate": kg_rate[i, :],
+                f"stack_{id}_curr_density": current_density[i, :],
+            }
+        )
+        stack_dfs.append(stack_df)
+
+    results_df = pd.concat([results_df, *stack_dfs], axis=1)
 
     return elec_sys, results_df
+
+
+def _run_electrolyzer_opt(modeling_options, power_signal):
+    # Tune to a desired system rating
+    options = calc_rated_system(modeling_options)
+
+    # Initialize system
+    elec_sys = Supervisor.from_dict(options["electrolyzer"])
+
+    # Define output variables
+    tot_kg = 0.0
+    max_curr_density = 0.0
+
+    # Run electrolyzer simulation
+    for i in range(len(power_signal)):
+        # TODO: replace with proper logging
+        # if (i % 1000) == 0:
+        #     print('Progress', i)
+        # print(i)
+        loop_H2, loop_h2_mfr, loop_power_left, curtailed = elec_sys.run_control(
+            power_signal[i]
+        )
+
+        tot_kg += loop_H2
+        new_curr = max([s.I / s.cell.cell_area for s in elec_sys.stacks])
+        max_curr_density = max(max_curr_density, new_curr)
+
+    return tot_kg, max_curr_density
+
+
+def run_electrolyzer(input_modeling, power_signal, optimize=False):
+    """
+    Runs an electrolyzer simulation based on a YAML configuration file and power
+    signal input.
+
+    Args:
+        input_modeling (`str` or `dict`): filepath specifying the YAML config
+            file, OR a dict representing a validated YAML config.
+        power_signal (`list`): An array representing power input
+        optimize (`bool`, optional): Whether the run will be based on an optimization.
+            For now, this entails tuning a system to a desired system rating, running
+            the simulation, and returning a simplified result for optimization runs.
+
+    Returns:
+        `Supervisor`: the instance used to run the simulation
+        `pandas.DataFrame`: a `DataFrame` representing the time series output
+            for the system, including values for each electrolyzer stack
+    """
+    err_msg = "Model input must be a str or dict object"
+    assert isinstance(
+        input_modeling,
+        (
+            str,
+            dict,
+        ),
+    ), err_msg
+
+    if isinstance(input_modeling, str):
+        # Parse/validate yaml configuration
+        modeling_options = val.load_modeling_yaml(input_modeling)
+    else:
+        modeling_options = input_modeling
+
+    if optimize:
+        return _run_electrolyzer_opt(modeling_options, power_signal)
+
+    return _run_electrolyzer_full(modeling_options, power_signal)

electrolyzer/inputs/modeling_schema.yaml

@@ -30,6 +30,14 @@ properties:
                 type: number
                 default: 1.0
                 description: simulation time step
+            initialize:
+                type: boolean
+                default: False
+                description: Determines whether the electrolyzer starts from an initial power (True), or from zero (False)
+            initial_power_kW:
+                type: number
+                default: 0.0
+                description: starting power for an initialized electrolyzer plant
 
             stack:
                 type: object
@@ -69,6 +77,10 @@ properties:
                 default: {}
                 description: Set control properties for electrolyzers
                 properties:
+                    system_rating_MW:
+                        type: number
+                        description: System rating
+                        unit: MW
                     n_stacks:
                         type: integer
                         default: 1