diff --git a/electrolyzer/components/classifiers/system_performance.py b/electrolyzer/components/classifiers/system_performance.py new file mode 100644 index 0000000..da1aad0 --- /dev/null +++ b/electrolyzer/components/classifiers/system_performance.py @@ -0,0 +1,42 @@ +import numpy as np +import openmdao.api as om + + +class SystemPerformance(om.ExplicitComponent): + """Connect clusters""" + + def initialize(self): + self.options.declare("n_clusters", types=(int, float), default=1.0) + + def setup(self): + self.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 self.vars_to_units.items(): + for i in range(0, int(self.options["n_clusters"])): + if ref_shape is None: + self.add_input(f"{v}_in_{i}", val=0.0, shape_by_conn=True, units=u) + ref_shape = f"{v}_in_{i}" + else: + self.add_input(f"{v}_in_{i}", val=0.0, copy_shape=ref_shape, units=u) + self.add_output(f"{v}_out", val=0.0, copy_shape=ref_shape, units=u) + + self.ref_shape = ref_shape + + def compute(self, inputs, outputs): + n_timesteps = len(inputs[self.ref_shape]) + + for v in self.vars_to_units.keys(): + var_cnt = np.zeros(n_timesteps) + for i in range(0, int(self.options["n_clusters"])): + var_cnt += inputs[f"{v}_in_{i}"] + outputs[f"{v}_out"] = var_cnt diff --git a/electrolyzer/components/cluster/simple_dynamics.py b/electrolyzer/components/cluster/simple_dynamics.py index 9db7958..6ade131 100644 --- a/electrolyzer/components/cluster/simple_dynamics.py +++ b/electrolyzer/components/cluster/simple_dynamics.py @@ -40,6 +40,5 @@ def compute(self, inputs, outputs): # But a separate output should be used to reflect show partial losses (warm-up delay) i_out = np.clip(inputs["I_in"], a_min=inputs["I_min"], a_max=inputs["I_max"]) - # TODO: add start-up delay - outputs["I_out"] = i_out + outputs["I_out"] = i_out * on_off_status outputs["on_off_status"] = on_off_status diff --git a/electrolyzer/control/openloop/simple_openloop_control.py b/electrolyzer/control/openloop/simple_openloop_control.py index 5c72859..c375658 100644 --- a/electrolyzer/control/openloop/simple_openloop_control.py +++ b/electrolyzer/control/openloop/simple_openloop_control.py @@ -14,10 +14,10 @@ def initialize(self): self.options.declare("plant_config", types=dict, default={}) self.options.declare("tech_config", types=dict, default={}) self.options.declare("n_clusters", types=int) - self.options.declare("control_variable", types=str, values=["power", "hydrogen"]) + self.options.declare("control_variable", values=["power", "hydrogen"]) def setup(self): - # self.n_timesteps = self.options["plant_config"]["simulation"]["n_timesteps"] + self.n_timesteps = self.options["plant_config"]["simulation"]["n_timesteps"] # self.dt = self.options["plant_config"]["simulation"]["dt"] # self.config = OLControlConfig.from_dict(self.options["tech_config"]["control_parameters"]) self.n_clusters = self.options["n_clusters"] @@ -25,12 +25,14 @@ def setup(self): if self.control_cmd == "power": # output_cmd_fmt = "power_cmd_{ci}" - self.add_input("P_command", val=0.0, shape_by_conn=True, units="kW") + # self.add_input("P_command", val=0.0, shape_by_conn=True, units="kW") + self.add_input("P_command", val=0.0, shape=self.n_timesteps, units="kW") for ci in range(self.n_clusters): self.add_output(f"P_command_{ci}", val=0.0, copy_shape="P_command", units="kW") else: # output_cmd_fmt = "hydrogen_cmd_{ci}" - self.add_input("H2_command", val=0.0, shape_by_conn=True, units="kg/h") + # self.add_input("H2_command", val=0.0, shape_by_conn=True, units="kg/h") + self.add_input("H2_command", val=0.0, shape=self.n_timesteps, units="kg/h") for ci in range(self.n_clusters): self.add_output(f"H2_command_{ci}", val=0.0, copy_shape="H2_command", units="kg/h") diff --git a/electrolyzer/core/bert.py b/electrolyzer/core/bert.py index 36d11f3..cdee455 100644 --- a/electrolyzer/core/bert.py +++ b/electrolyzer/core/bert.py @@ -12,9 +12,7 @@ from electrolyzer.connectors.degradation_combiner import CombineDegradation 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 +from electrolyzer.components.classifiers.system_performance import SystemPerformance class State(IntEnum): @@ -25,23 +23,31 @@ class State(IntEnum): class BERT: - def __init__(self, config_input, make_n2=True): + def __init__(self, config_input, make_n2=True, as_problem=True): self.create_n2 = make_n2 self.supported_models = supported_models.copy() # read in config file; it's a yaml dict that looks like this: self.load_config(config_input) - self.prob = om.Problem(reports=False) - self.model = self.prob.model - plant_group = om.Group() - # Create the plant model group and add components - self.plant = self.model.add_subsystem("plant", plant_group, promotes=["*"]) + if as_problem: + self.prob = om.Problem(reports=False) + self.model = self.prob.model + plant_group = om.Group() + + # Create the plant model group and add components + self.plant = self.model.add_subsystem("plant", plant_group, promotes=["*"]) + else: + self.plant = om.Group() self.create_controller() self.create_components() + self.create_performance_aggregator() + + self.connect_system() - self.create_recorder(self.prob) + if as_problem: + self.create_recorder(self.prob) self.state = State.INITIALIZED @@ -98,31 +104,11 @@ def run(self): def post_process(self): pass - def create_cluster_components(self): - pass - - def create_controller(self): - controller = self.create_controller_component() - self.plant.add_subsystem("controller", controller) - - def create_components(self): - # - + def create_cluster_group(self): # Get the design parameters of the cell cell_design_params = get_cell_params_for_model(self.config["cell"].get("model", None)) - # Step 1: Create cluster groups - clusters = [] - cluster_i = 0 - - # NOTE: cell design params should only be promoted if all the clusters are identical - if self.identical_cells: - cluster_group = self.plant.add_subsystem( - f"Cluster{cluster_i}", om.Group(), promotes=cell_design_params - ) - else: - cluster_group = self.plant.add_subsystem(f"Cluster{cluster_i}", om.Group()) - clusters.append(cluster_group) + cluster_group = om.Group() # Step 2: Create controller cluster connector components pre_translator = self.create_controller_cluster_connector(cell_design_params) @@ -156,16 +142,76 @@ def create_components(self): f"converter.ref_cell.{var}_cell_out", f"classifier.cell_classifier.{var}_in" ) + return cluster_group + + def create_controller(self): + controller = self.create_controller_component() + self.plant.add_subsystem("controller", controller) + + def create_performance_aggregator(self): + ts_perf_mod = SystemPerformance(n_clusters=self.n_clusters) + self.plant.add_subsystem("system_timeseries", ts_perf_mod) + + perf_mod = SystemPerformance(n_clusters=self.n_clusters) + self.plant.add_subsystem("system_ub", perf_mod) + + def create_components(self): + # + + # Get the design parameters of the cell + cell_design_params = get_cell_params_for_model(self.config["cell"].get("model", None)) + stack_design_params = [*cell_design_params, "n_stacks", "n_cells"] + # Step 1: Create cluster groups + clusters = [] + # cluster_i = 0 + for cluster_i in range(self.n_clusters): + cluster_comp = self.create_cluster_group() + + # NOTE: cell design params should only be promoted if all the clusters are identical + if self.identical_cells: + cluster_group = self.plant.add_subsystem( + f"Cluster{cluster_i}", cluster_comp, promotes=stack_design_params + ) + else: + cluster_group = self.plant.add_subsystem(f"Cluster{cluster_i}", cluster_comp) + clusters.append(cluster_group) + # Connect controller to cluster - self.plant.connect( - "controller.P_command", f"Cluster{cluster_i}.translator.cluster_to_stack.P_in" - ) + # self.plant.connect( + # f"controller.{self.control_passed_var}_command_{cluster_i}", + # f"Cluster{cluster_i}.translator.cluster_to_stack.{self.control_passed_var}_in" + # ) # cluster_group.connect("converter.I_ref_points", "classifier." # cluster_group.connect("converter.ref_cell.") self.clusters = clusters + def connect_system(self): + # Connect controller to cluster + + for cluster_i in range(0, self.n_clusters, 1): + # Connect controller to cluster + self.plant.connect( + f"controller.{self.control_passed_var}_command_{cluster_i}", + f"Cluster{cluster_i}.translator.cluster_to_stack.{self.control_passed_var}_in", + ) + + for cluster_i in range(0, self.n_clusters, 1): + # connect the clusters to a system performance component + # connect the classifier component and the simulation component + for var in ["P", "H2", "O2", "V"]: + self.plant.connect( + # part of scale_stack_to_cluster + f"Cluster{cluster_i}.simulation.Cluster_{var}", + f"system_timeseries.{var}_in_{cluster_i}", + ) + self.plant.connect( + # part of classifier.stack_to_cluster_ub + f"Cluster{cluster_i}.classifier.{var}_max", + f"system_ub.{var}_in_{cluster_i}", + ) + def create_cluster_simulation_block(self, cell_design_params): simulation = om.Group() @@ -416,11 +462,19 @@ def create_controller_component(self): n_timesteps = int(self.plant_config["simulation"]["n_timesteps"]) if "control_model" not in self.system_config: ivc_comp = om.IndepVarComp( - name=f"{self.control_passed_var}_command", val=np.full(n_timesteps, 40.0), units="W" + name=f"{self.control_passed_var}_command_0", + val=np.full(n_timesteps, 40.0), + units="W", ) + if self.n_clusters > 1: + msg = ( + "Cannot run multiple clusters without a control model. " + "Please specify a control model" + ) + raise NotImplementedError(msg) return ivc_comp controller_name = self.system_config["control_model"] - controller_model = self.supported_models(controller_name) + controller_model = self.supported_models.get(controller_name) controller = controller_model( plant_config=self.plant_config, tech_config=self.system_config, diff --git a/electrolyzer/test/test_om_examples.py b/electrolyzer/test/test_om_examples.py index 7dd135f..cfa7726 100644 --- a/electrolyzer/test/test_om_examples.py +++ b/electrolyzer/test/test_om_examples.py @@ -6,6 +6,7 @@ from electrolyzer import BERT_EXAMPLE_DIR from electrolyzer.core.bert import BERT +from electrolyzer.core.file_utils import load_yaml def test_example_00_no_controller(subtests): @@ -13,16 +14,18 @@ def test_example_00_no_controller(subtests): os.chdir(example_fpath) config_fpath = example_fpath / "bert_config.yaml" - bert = BERT(config_fpath, make_n2=False) + config = load_yaml(config_fpath) + config["simulation"]["n_timesteps"] = 20 + bert = BERT(config, make_n2=False) bert.run() - scale_fac = bert.model.get_val("Cluster0.n_stacks", units="unitless") * bert.model.get_val( - "Cluster0.n_cells", units="unitless" + scale_fac = bert.model.get_val("n_stacks", units="unitless") * bert.model.get_val( + "n_cells", units="unitless" ) p_cell_ref = bert.model.get_val("Cluster0.converter.ref_cell.P_cell_out", units="W") p_system_ref = p_cell_ref * scale_fac - bert.model.set_val("controller.P_command", p_system_ref, units="W") + bert.model.set_val("controller.P_command_0", p_system_ref, units="W") bert.run() i_estimated = bert.model.get_val("Cluster0.translator.command_to_current.I_command", units="A") i_actual = bert.model.get_val("Cluster0.converter.I_ref_points", units="A") @@ -30,10 +33,10 @@ def test_example_00_no_controller(subtests): i_error = i_estimated - i_actual with subtests.test("100 cells per stack"): - assert pytest.approx(100.0, rel=1e-6) == bert.model.get_val("Cluster0.n_cells")[0] + assert pytest.approx(100.0, rel=1e-6) == bert.model.get_val("n_cells")[0] with subtests.test("1 stack per cluster"): - assert pytest.approx(1.0, rel=1e-6) == bert.model.get_val("Cluster0.n_stacks")[0] + assert pytest.approx(1.0, rel=1e-6) == bert.model.get_val("n_stacks")[0] with subtests.test("I-V Curve fit error is less than 0.105 A"): assert np.all(np.abs(i_error) < 0.105) @@ -117,14 +120,14 @@ def test_example_00_no_controller(subtests): ) assert ( pytest.approx( - cell_rated_power * bert.model.get_val("Cluster0.n_cells", units="unitless"), + cell_rated_power * bert.model.get_val("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") + ) == stack_rated_power * bert.model.get_val("n_stacks", units="unitless") with subtests.test("Rated conversion efficiency"): assert pytest.approx(60.84498639, rel=1e-6) == bert.model.get_val( @@ -187,3 +190,86 @@ def test_example_00_no_controller(subtests): P_cell_deg = I_deg * (V_cell_deg + V_cell_bol) assert np.allclose(P_cell_bol, P_cell_deg) assert np.allclose(P_cell_bol * scale_fac, P_cell_deg * scale_fac) + + +def test_example_00_with_controller(subtests): + example_fpath = BERT_EXAMPLE_DIR / "example_00_refactor" + os.chdir(example_fpath) + + n_clusters = 2 + + config_fpath = example_fpath / "bert_config.yaml" + config = load_yaml(config_fpath) + config["simulation"]["n_timesteps"] = 20 + config["system"]["n_clusters"] = n_clusters + config["system"]["control_model"] = "OLBasicSplit" + bert = BERT(config, make_n2=False) + bert.run() + + scale_fac = bert.model.get_val("n_stacks", units="unitless") * bert.model.get_val( + "n_cells", units="unitless" + ) + + p_cell_ref = bert.model.get_val("Cluster0.converter.ref_cell.P_cell_out", units="W") + p_system_ref = p_cell_ref * scale_fac * n_clusters + bert.model.set_val("controller.P_command", p_system_ref, units="W") + bert.run() + + with subtests.test("Initial reference rated power"): + assert ( + pytest.approx(4467.1560, rel=1e-6) + == bert.model.get_val("Cluster0.converter.ref_cell.P_cell_out", units="W")[-1] + ) + + coeff_initial = copy.deepcopy( + bert.model.get_val("Cluster0.translator.command_to_current.curve_coeffs", units="A/W") + ) + + with subtests.test("Initial curve coefficients"): + expected_initial_coeff = np.array( + [7.08472908e-10, -1.70727901e-05, 4.78002528e-01, 2.34225327e00, -1.42414827e01] + ) + assert pytest.approx(expected_initial_coeff, rel=1e-6, abs=1e-8) == coeff_initial + with subtests.test("Initial reference rated voltage"): + assert ( + pytest.approx(2.233578003273652, rel=1e-6) + == bert.model.get_val("Cluster0.converter.ref_cell.V_cell_out", units="V")[-1] + ) + # TODO: also test Cluster1 + + with subtests.test("Cluster0 on/off status"): + assert ( + 19.0 + == bert.model.get_val( + "Cluster0.simulation.dynamics.on_off_status", units="unitless" + ).sum() + ) + + with subtests.test("Cluster1 on/off status"): + assert ( + 19.0 + == bert.model.get_val( + "Cluster1.simulation.dynamics.on_off_status", units="unitless" + ).sum() + ) + + with subtests.test("Cluster H2 production"): + assert ( + pytest.approx(81.64025077453496, rel=1e-6) + == bert.model.get_val("Cluster0.simulation.Cluster_H2", units="kg/h").sum() + ) + assert ( + pytest.approx(81.64025077453496, rel=1e-6) + == bert.model.get_val("Cluster1.simulation.Cluster_H2", units="kg/h").sum() + ) + + with subtests.test("System H2 Production"): + assert ( + pytest.approx(163.28050154906992, rel=1e-6) + == bert.model.get_val("system_timeseries.H2_out", units="kg/h").sum() + ) + + with subtests.test("System Rated H2 Production"): + assert pytest.approx(14.982832485661488, rel=1e-6) == bert.model.get_val( + "system_ub.H2_out", units="kg/h" + ) diff --git a/examples/example_00_refactor/bert_config.yaml b/examples/example_00_refactor/bert_config.yaml index 472aeba..bc92f93 100644 --- a/examples/example_00_refactor/bert_config.yaml +++ b/examples/example_00_refactor/bert_config.yaml @@ -1,6 +1,6 @@ simulation: dt: 3600 # seconds - n_timesteps: 20 # number of timesteps, UNUSED + n_timesteps: 8760 # number of timesteps # electrolyzer: system: n_clusters: 1