diff --git a/electrolyzer/components/classifiers/__init__.py b/electrolyzer/components/classifiers/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/electrolyzer/components/classifiers/cell_classifier.py b/electrolyzer/components/classifiers/cell_classifier.py new file mode 100644 index 0000000..a5e1f04 --- /dev/null +++ b/electrolyzer/components/classifiers/cell_classifier.py @@ -0,0 +1,56 @@ +import numpy as np +import openmdao.api as om + + +class CellClassification(om.ExplicitComponent): + def initialize(self): + self.options.declare("tech_config", types=dict, default={}) + self.options.declare("plant_config", types=dict, default={}) + + def setup(self): + self.add_input("I_max", val=0.0, shape=1, units="A") + self.add_input("I_min", val=0.0, shape=1, units="A") + self.add_input("I_ref_points", val=0.0, shape_by_conn=True, units="A") + + self.vars_to_units = { + "J": "A/(cm**2)", + "P": "kW", + "H2": "kg/h", + "O2": "kg/h", + # "H2O": "kg/s", + "V": "V", + } + + ref_shape = "I_ref_points" + for v, u in self.vars_to_units.items(): + self.add_input(f"{v}_in", val=0.0, copy_shape=ref_shape, units=u) + self.add_output(f"{v}_min", val=0.0, shape=1, units=u) + self.add_output(f"{v}_max", val=0.0, shape=1, units=u) + + # energy bounds + self.add_output("efficiency_min", val=0.0, shape=1, units="kW*h/kg") + self.add_output("efficiency_max", val=0.0, shape=1, units="kW*h/kg") + # Should output: + # - rated cell voltage + # - rated current density + # - rated power consumption + # - rated h2 production rate + # - rated efficiency + # - rated o2 production rate + # - rated water consumption rate + + # self.add_input("J_max", val=self.config.J_max, shape=1, units="A/(cm**2)") + # self.add_input("J_min", val=self.config.J_max, shape=1, units="A/(cm**2)") + + def compute(self, inputs, outputs): + idx_ref_min = np.argwhere(inputs["I_ref_points"] <= inputs["I_min"]).flatten()[-1] + idx_ref_max = np.argwhere(inputs["I_ref_points"] >= inputs["I_max"]).flatten()[0] + + for v in list(self.vars_to_units.keys()): + outputs[f"{v}_min"] = inputs[f"{v}_in"][idx_ref_min] + outputs[f"{v}_max"] = inputs[f"{v}_in"][idx_ref_max] + + # kWh/kg + efficiency = inputs["P_in"] / inputs["H2_in"] + outputs["efficiency_min"] = efficiency[idx_ref_min] + outputs["efficiency_max"] = efficiency[idx_ref_max] diff --git a/electrolyzer/connectors/series_scalar.py b/electrolyzer/connectors/series_scalar.py new file mode 100644 index 0000000..e47f1d1 --- /dev/null +++ b/electrolyzer/connectors/series_scalar.py @@ -0,0 +1,64 @@ +import openmdao.api as om + + +class GenericSeriesConverter(om.ExplicitComponent): + def initialize(self): + self.options.declare("scaling_component", types=str) + self.options.declare("n_comps", types=(int, float), default=1.0) + + def setup(self): + self.add_input( + f"n_{self.options['scaling_component']}", + val=self.options["n_comps"], + shape=1, + units="unitless", + ) + vars_to_units = { + "J": "A/(cm**2)", + "I": "A", + "P": "W", + "H2": "kg/s", + "O2": "kg/s", + # "H2O": "kg/s", + "V": "V", + # "V_deg": "V" + } + + ref_shape = None + for v, u in vars_to_units.items(): + if ref_shape is None: + self.add_input(f"{v}_in", val=0.0, shape_by_conn=True, units=u) + self.add_output(f"{v}_out", val=0.0, copy_shape=f"{v}_in", units=u) + ref_shape = f"{v}_in" + else: + self.add_input(f"{v}_in", val=0.0, copy_shape=ref_shape, units=u) + self.add_output(f"{v}_out", val=0.0, copy_shape=ref_shape, units=u) + + +class SplitAcrossSerialComponents(GenericSeriesConverter): + """Scale power down""" + + def compute(self, inputs, outputs): + for o_name in outputs.keys(): + in_name = o_name.replace("_out", "_in") + + if o_name == "I_out" or o_name == "J_out": + outputs[o_name] = inputs[in_name] + else: + outputs[o_name] = inputs[in_name] / inputs[f"n_{self.options['scaling_component']}"] + + +class CombineSerialComponents(GenericSeriesConverter): + """Scale power down""" + + def setup(self): + super().setup() + + def compute(self, inputs, outputs): + for o_name in outputs.keys(): + in_name = o_name.replace("_out", "_in") + + if o_name == "I_out" or o_name == "J_out": + outputs[o_name] = inputs[in_name] + else: + outputs[o_name] = inputs[in_name] * inputs[f"n_{self.options['scaling_component']}"] diff --git a/electrolyzer/core/bert.py b/electrolyzer/core/bert.py index fb64f7a..01856eb 100644 --- a/electrolyzer/core/bert.py +++ b/electrolyzer/core/bert.py @@ -4,9 +4,16 @@ import numpy as np import openmdao.api as om -from electrolyzer.core.file_utils import load_yaml +from electrolyzer.core.file_utils import load_yaml, make_unique_case_name from electrolyzer.core.supported_models import supported_models +from electrolyzer.connectors.series_scalar import ( + CombineSerialComponents, # , SplitAcrossSerialComponents +) from electrolyzer.components.cell.cell_design_params import get_cell_params_for_model +from electrolyzer.components.classifiers.cell_classifier import CellClassification + + +# from electrolyzer.components.classifiers.system_performance import SystemPerformance class State(IntEnum): @@ -33,6 +40,8 @@ def __init__(self, config_input, make_n2=True): self.create_controller() self.create_components() + self.create_recorder(self.prob) + self.state = State.INITIALIZED def load_config(self, config_input): @@ -116,6 +125,7 @@ def create_components(self): # Step 2: Create controller cluster connector components pre_translator = self.create_controller_cluster_connector(cell_design_params) + cluster_classifier = self.create_cluster_classification_component() # Translator has scale down + power to current conversion translator = self.create_controller_translator() # Step 3: Create the simulate block of a cluster @@ -123,6 +133,9 @@ def create_components(self): promotion_vars = [*cell_design_params, "I_min", "I_max"] cluster_group.add_subsystem("converter", pre_translator, promotes=promotion_vars) + cluster_group.add_subsystem( + "classifier", cluster_classifier, promotes=["I_min", "I_max", "n_cells", "n_stacks"] + ) cluster_group.add_subsystem("translator", translator, promotes=["n_stacks", "n_cells"]) cluster_group.add_subsystem("simulation", simulator, promotes=promotion_vars) @@ -132,10 +145,23 @@ def create_components(self): "converter.p2i.curve_coeffs", "translator.command_to_current.curve_coeffs" ) cluster_group.connect("translator.command_to_current.I_command", "simulation.dynamics.I_in") + + # connect converter group stuff to classification block + cluster_group.connect("converter.I_ref_points", "classifier.I_ref_points") + cluster_group.connect("converter.ref_cell.J_out", "classifier.cell_classifier.J_in") + for var in ["P", "H2", "O2", "V"]: + cluster_group.connect( + f"converter.ref_cell.{var}_cell_out", f"classifier.cell_classifier.{var}_in" + ) + + # Connect controller to cluster self.plant.connect( "controller.P_command", f"Cluster{cluster_i}.translator.cluster_to_stack.P_in" ) + # cluster_group.connect("converter.I_ref_points", "classifier." + # cluster_group.connect("converter.ref_cell.") + self.clusters = clusters def create_cluster_simulation_block(self, cell_design_params): @@ -198,8 +224,22 @@ def create_controller_translator(self): return translator def create_controller_cluster_connector(self, cell_design_params): + """Group containing the: + + 1. Bounds and reference point component (min/max bounds) + 2. Reference cell component + 3. Curve coefficient component + + Args: + cell_design_params (list[str]): cell design parameters to promote + + Returns: + om.Group: pre-simulation group + """ pre_converter_grp = om.Group() bounds_comp = self.create_bounds_component() + + # 1. Operational bounds and reference point component pre_converter_grp.add_subsystem( "IJ_ref", bounds_comp, @@ -207,19 +247,89 @@ def create_controller_cluster_connector(self, cell_design_params): promotes_outputs=["I_ref_points", "I_min", "I_max"], ) + # 2. Reference cell model cell = self.create_cell_model() pre_converter_grp.add_subsystem("ref_cell", cell, promotes_inputs=cell_design_params) + # 3. Curve coefficient component coeff_comp = self.create_component("control_command_converter", model_key="coeff_model") pre_converter_grp.add_subsystem("p2i", coeff_comp, promotes_inputs=["I_ref_points"]) + # Connect components + # Connect the reference points to the cell pre_converter_grp.connect("I_ref_points", "ref_cell.I_in") # Connect the power output from the cell to the power to current thing pre_converter_grp.connect("ref_cell.P_cell_out", "p2i.P_ref_points") - return pre_converter_grp + def create_cluster_classification_component(self): + # 4. Add the cell classification component to the system + # The cell classification component inputs of J_in, P_in, H2_in, O2_in, V_in + # and outputs min and max values of each input, plus min/max efficiency values + # Cell model outputs P_cell_out, J_out, H2_cell_out, O2_cell_out, V_cell_out + classifier_group = om.Group() + + cell_classifier = CellClassification(tech_config={}, plant_config=self.plant_config) + classifier_group.add_subsystem( + "cell_classifier", + cell_classifier, + promotes_inputs=["I_ref_points", "I_min", "I_max"], + promotes_outputs=["efficiency_min", "efficiency_max"], + ) + + cell_scale_up_lb = CombineSerialComponents( + scaling_component="cells", n_comps=self.config["stack"]["n_cells"] + ) + stack_scale_up_lb = CombineSerialComponents( + scaling_component="stacks", n_comps=self.config["cluster"]["n_stacks"] + ) + + cell_scale_up_ub = CombineSerialComponents( + scaling_component="cells", n_comps=self.config["stack"]["n_cells"] + ) + stack_scale_up_ub = CombineSerialComponents( + scaling_component="stacks", n_comps=self.config["cluster"]["n_stacks"] + ) + + bounds_base_vars = ["J", "P", "H2", "O2", "V"] + + promoted_outputs_lb = [(f"{v}_out", f"{v}_min") for v in bounds_base_vars] + promoted_outputs_ub = [(f"{v}_out", f"{v}_max") for v in bounds_base_vars] + + # # lower bounds + classifier_group.add_subsystem( + "cell_to_stack_lb", cell_scale_up_lb, promotes_inputs=["n_cells"] + ) + classifier_group.add_subsystem( + "stack_to_cluster_lb", + stack_scale_up_lb, + promotes_inputs=["n_stacks"], + promotes_outputs=promoted_outputs_lb, + ) + # # upper bounds + classifier_group.add_subsystem( + "cell_to_stack_ub", cell_scale_up_ub, promotes_inputs=["n_cells"] + ) + classifier_group.add_subsystem( + "stack_to_cluster_ub", + stack_scale_up_ub, + promotes_inputs=["n_stacks"], + promotes_outputs=promoted_outputs_ub, + ) + # Cluster0.classifier.cell_to_stack_lb.I_in + # classifier_group.connect("I_min", "cell_to_stack_lb.I_in") + for var in bounds_base_vars: + # # scale-up lower bounds + classifier_group.connect(f"cell_classifier.{var}_min", f"cell_to_stack_lb.{var}_in") + classifier_group.connect(f"cell_to_stack_lb.{var}_out", f"stack_to_cluster_lb.{var}_in") + + # # scale-up upper bounds + classifier_group.connect(f"cell_classifier.{var}_max", f"cell_to_stack_ub.{var}_in") + classifier_group.connect(f"cell_to_stack_ub.{var}_out", f"stack_to_cluster_ub.{var}_in") + + return classifier_group + def create_cell_model(self): cell_config = self.config["cell"] if (cell_model_name := cell_config.get("model", None)) is not None: @@ -275,3 +385,97 @@ def create_controller_component(self): control_variable=self.control_var, ) return controller + + def create_recorder(self, opt_prob): + # TODO: put this into pose_optimization one day + + if "recorder" not in self.config: + return None + + folder_output = self.config.get("folder_output", Path.cwd()) + + recorder_options = ["record_inputs", "record_outputs", "record_residuals"] + if self.config["recorder"].get("flag", False): + # Check that the output folder exists and create it if needed + if not Path(folder_output).exists(): + Path.mkdir(folder_output, parents=True, exist_ok=True) + + if self.config["recorder"].get("flag", False): + # Check that the output folder exists and create it if needed + if not Path(folder_output).exists(): + Path.mkdir(folder_output, parents=True, exist_ok=True) + + overwrite_recorder = self.config["recorder"].get("overwrite_recorder", False) + recorder_path = Path(folder_output) / self.config["recorder"]["file"] + + if not overwrite_recorder: + # make a unique filename with the same base as self.config["recorder"]["file"] + # separate out the filename without the extension + file_base = self.config["recorder"]["file"].split(".sql")[0] + + recorder_fname = make_unique_case_name( + Path(folder_output), f"{file_base}.sql", ".sql" + ) + recorder_path = Path(folder_output) / recorder_fname + + recorder_attachment = ( + self.config["recorder"].get("recorder_attachment", "driver").lower() + ) + allowed_attachments = ["driver", "model"] + if recorder_attachment not in allowed_attachments: + msg = ( + f"Invalid recorder attachment '{recorder_attachment}'. " + f"Currently supported options are {allowed_attachments}. " + "We recommend using 'driver' if running an optimization " + "or parameter sweep in parallel." + ) + raise ValueError(msg) + + # Create recorder + recorder = om.SqliteRecorder(recorder_path) + + if recorder_attachment == "model": + # add the recorder to the model + recorder_options += ["options_excludes"] + + opt_prob.model.add_recorder(recorder) + + for recorder_opt in recorder_options: + if recorder_opt in self.config["recorder"]: + opt_prob.model.recording_options[recorder_opt] = self.config[ + "recorder" + ].get(recorder_opt) + + opt_prob.model.recording_options["includes"] = self.config["recorder"].get( + "includes", ["*"] + ) + # opt_prob.model.recording_options["excludes"] = self.config["recorder"].get( + # "excludes", ["*resource_data"] + # ) + return recorder_path + + if recorder_attachment == "driver": + recorder_options += [ + "record_constraints", + "record_derivative", + "record_desvars", + "record_objectives", + ] + # add the recorder to the driver + opt_prob.driver.add_recorder(recorder) + + for recorder_opt in recorder_options: + if recorder_opt in self.config["recorder"]: + opt_prob.driver.recording_options[recorder_opt] = self.config[ + "recorder" + ].get(recorder_opt) + + opt_prob.driver.recording_options["includes"] = self.config["recorder"].get( + "includes", ["*"] + ) + # opt_prob.driver.recording_options["excludes"] = self.config["recorder"].get( + # "excludes", ["*resource_data"] + # ) + return recorder_path + + return None diff --git a/electrolyzer/core/file_utils.py b/electrolyzer/core/file_utils.py index 26dab44..cc4a8bc 100644 --- a/electrolyzer/core/file_utils.py +++ b/electrolyzer/core/file_utils.py @@ -1,3 +1,4 @@ +import re from pathlib import Path import yaml @@ -8,3 +9,43 @@ def load_yaml(filename, loader=yaml.SafeLoader) -> dict: return filename # filename already yaml dict with Path.open(filename) as fid: return yaml.load(fid, loader) + + +def make_unique_case_name(folder, proposed_fname, fext): + """Generate a filename that does not already exist in a user-defined folder. + + Args: + folder (str | Path): directory that a file is expected to be created in. + proposed_fname (str): filename (with extension) to check for existence and + to use as the base file description of a new an unique file name. + fext (str): file extension, such as ".csv", ".sql", ".yaml", etc. + + Returns: + str: unique filename that does not yet exist in folder. + """ + if "." not in fext: + fext = f".{fext}" + + # if file(s) exist with the same base name, make a new unique filename + file_base = proposed_fname.split(fext)[0] + existing_files = [f for f in Path(folder).glob(f"**/*{fext}") if file_base in f.name] + if len(existing_files) == 0: + return proposed_fname + + # get past numbers that were used to make unique files by matching + # filenames against the file base name followed by a number + past_numbers = [ + int(re.findall(f"{file_base}[0-9]+", str(fname))[0].split(file_base)[-1]) + for fname in existing_files + if len(re.findall(f"{file_base}[0-9]+", str(fname))) > 0 + ] + + if len(past_numbers) > 0: + # if multiple files have the same basename followed by a number, + # take the maximum unique number and add one + unique_number = int(max(past_numbers) + 1) + return f"{file_base}{unique_number}{fext}" + else: + # if no files have the same basename followed by a number, + # but do have the same basename, then add a zero to the file basename + return f"{file_base}0{fext}" diff --git a/electrolyzer/test/test_om_examples.py b/electrolyzer/test/test_om_examples.py index 3bc5ddc..0c5b66d 100644 --- a/electrolyzer/test/test_om_examples.py +++ b/electrolyzer/test/test_om_examples.py @@ -99,3 +99,45 @@ def test_example_00_no_controller(subtests): [8.30333109e-10, -1.92246906e-05, 4.75331901e-01, 2.52511351e00, -1.62209719e01] ) assert pytest.approx(expected_coeff, rel=1e-6, abs=1e-8) == coeff_new + + with subtests.test("Cluster min voltage"): + assert ( + pytest.approx( + bert.model.get_val("Cluster0.classifier.cell_classifier.V_max", units="V"), rel=1e-6 + ) + == bert.model.get_val("Cluster0.classifier.V_max", units="V") / scale_fac + ) + + with subtests.test("Cell/stack/cluster max power"): + cell_rated_power = bert.model.get_val( + "Cluster0.classifier.cell_to_stack_ub.P_in", units="kW" + ) + stack_rated_power = bert.model.get_val( + "Cluster0.classifier.stack_to_cluster_ub.P_in", units="kW" + ) + assert ( + pytest.approx( + cell_rated_power * bert.model.get_val("Cluster0.n_cells", units="unitless"), + rel=1e-6, + ) + == stack_rated_power + ) + assert pytest.approx( + bert.model.get_val("Cluster0.classifier.P_max", units="kW"), rel=1e-6 + ) == stack_rated_power * bert.model.get_val("Cluster0.n_stacks", units="unitless") + + with subtests.test("Rated conversion efficiency"): + assert pytest.approx(60.84498639, rel=1e-6) == bert.model.get_val( + "Cluster0.classifier.efficiency_max", units="kW*h/kg" + ) + + with subtests.test("Min conversion efficiency"): + assert pytest.approx(48.23406508, rel=1e-6) == bert.model.get_val( + "Cluster0.classifier.efficiency_min", units="kW*h/kg" + ) + + with subtests.test("H2 rated production"): + assert ( + pytest.approx(7.491416242830744, rel=1e-6) + == bert.model.get_val("Cluster0.classifier.H2_max", units="kg/h")[0] + ) diff --git a/examples/example_00_refactor/run.py b/examples/example_00_refactor/run.py index 57e8df0..79d9f11 100644 --- a/examples/example_00_refactor/run.py +++ b/examples/example_00_refactor/run.py @@ -2,9 +2,24 @@ from pathlib import Path from electrolyzer.core.bert import BERT +from electrolyzer.core.file_utils import load_yaml os.chdir(Path(__file__).parent) config_fpath = Path(__file__).parent / "bert_config.yaml" bert = BERT(config_fpath) bert.run() + + +# Run with recorder +config = load_yaml(config_fpath) +config["folder_output"] = Path(__file__).parent / "outputs" +config["recorder"] = { + "flag": True, + "file": "case.sql", + "overwrite_recorder": True, + "recorder_attachment": "model", + "includes": ["*"], +} +bert = BERT(config) +bert.run()