Allow input parameters in ESOH model

CHANGELOG.md

@@ -18,6 +18,7 @@
 - Updated `_steps_util.py` to throw a specific exception when drive cycle starts at t>0 ([#3756](https://github.com/pybamm-team/PyBaMM/pull/3756))
 - Updated `plot_voltage_components.py` to support both `Simulation` and `Solution` objects. Added new methods in both `Simulation` and `Solution` classes for allow the syntax `simulation.plot_voltage_components` and `solution.plot_voltage_components`. Updated `test_plot_voltage_components.py` to reflect these changes ([#3723](https://github.com/pybamm-team/PyBaMM/pull/3723)).
 - The SEI thickness decreased at some intervals when the 'electron-migration limited' model was used. It has been corrected ([#3622](https://github.com/pybamm-team/PyBaMM/pull/3622))
+- Allow input parameters in ESOH model ([#3921](https://github.com/pybamm-team/PyBaMM/pull/3921))
 
 ## Optimizations
 

pybamm/models/full_battery_models/lithium_ion/electrode_soh.py

@@ -380,7 +380,8 @@ def solve(self, inputs):
         # Calculate theoretical energy
         # TODO: energy calc for MSMR
         if self.options["open-circuit potential"] != "MSMR":
-            energy = self.theoretical_energy_integral(sol_dict)
+            energy_inputs = {**sol_dict, **inputs}
+            energy = self.theoretical_energy_integral(energy_inputs)
             sol_dict.update({"Maximum theoretical energy [W.h]": energy})
         return sol_dict
 
@@ -604,12 +605,10 @@ def _check_esoh_feasible(self, inputs):
         else:
             # address numpy 1.25 deprecation warning: array should have ndim=0
             # before conversion
-            V_lower_bound = float(
-                self.OCV_function.evaluate(inputs={"x": x0_min, "y": y0_max}).item()
-            )
-            V_upper_bound = float(
-                self.OCV_function.evaluate(inputs={"x": x100_max, "y": y100_min}).item()
-            )
+            inputs.update({"x": x0_min, "y": y0_max}, check_already_exists=False)
+            V_lower_bound = float(self.OCV_function.evaluate(inputs=inputs).item())
+            inputs.update({"x": x100_max, "y": y100_min})
+            V_upper_bound = float(self.OCV_function.evaluate(inputs=inputs).item())
 
         # Check that the min and max achievable voltages span wider than the desired
         # voltage range
@@ -632,7 +631,7 @@ def _check_esoh_feasible(self, inputs):
                 )
             )
 
-    def get_initial_stoichiometries(self, initial_value, tol=1e-6):
+    def get_initial_stoichiometries(self, initial_value, tol=1e-6, inputs=None):
         """
         Calculate initial stoichiometries to start off the simulation at a particular
         state of charge, given voltage limits, open-circuit potentials, etc defined by
@@ -657,7 +656,7 @@ def get_initial_stoichiometries(self, initial_value, tol=1e-6):
         """
         parameter_values = self.parameter_values
         param = self.param
-        x_0, x_100, y_100, y_0 = self.get_min_max_stoichiometries()
+        x_0, x_100, y_100, y_0 = self.get_min_max_stoichiometries(inputs=inputs)
 
         if isinstance(initial_value, str) and initial_value.endswith("V"):
             V_init = float(initial_value[:-1])
@@ -714,7 +713,7 @@ def get_initial_stoichiometries(self, initial_value, tol=1e-6):
 
         return x, y
 
-    def get_min_max_stoichiometries(self):
+    def get_min_max_stoichiometries(self, inputs=None):
         """
         Calculate min/max stoichiometries
         given voltage limits, open-circuit potentials, etc defined by parameter_values
@@ -724,18 +723,23 @@ def get_min_max_stoichiometries(self):
         x_0, x_100, y_100, y_0
             The min/max stoichiometries
         """
+        inputs = inputs or {}
         parameter_values = self.parameter_values
         param = self.param
 
-        Q_n = parameter_values.evaluate(param.n.Q_init)
-        Q_p = parameter_values.evaluate(param.p.Q_init)
+        Q_n = parameter_values.evaluate(param.n.Q_init, inputs=inputs)
+        Q_p = parameter_values.evaluate(param.p.Q_init, inputs=inputs)
 
         if self.known_value == "cyclable lithium capacity":
-            Q_Li = parameter_values.evaluate(param.Q_Li_particles_init)
-            inputs = {"Q_n": Q_n, "Q_p": Q_p, "Q_Li": Q_Li}
+            Q_Li = parameter_values.evaluate(param.Q_Li_particles_init, inputs=inputs)
+            inputs.update(
+                {"Q_n": Q_n, "Q_p": Q_p, "Q_Li": Q_Li}, check_already_exists=False
+            )
         elif self.known_value == "cell capacity":
-            Q = parameter_values.evaluate(param.Q / param.n_electrodes_parallel)
-            inputs = {"Q_n": Q_n, "Q_p": Q_p, "Q": Q}
+            Q = parameter_values.evaluate(
+                param.Q / param.n_electrodes_parallel, inputs=inputs
+            )
+            inputs.update({"Q_n": Q_n, "Q_p": Q_p, "Q": Q}, check_already_exists=False)
         # Solve the model and check outputs
         sol = self.solve(inputs)
         return [sol["x_0"], sol["x_100"], sol["y_100"], sol["y_0"]]
@@ -814,7 +818,7 @@ def theoretical_energy_integral(self, inputs, points=1000):
         param = self.param
         T = param.T_amb_av(0)
         Vs = self.parameter_values.evaluate(
-            param.p.prim.U(y_vals, T) - param.n.prim.U(x_vals, T)
+            param.p.prim.U(y_vals, T) - param.n.prim.U(x_vals, T), inputs=inputs
         ).flatten()
         # Calculate dQ
         Q = Q_p * (y_0 - y_100)
@@ -831,6 +835,7 @@ def get_initial_stoichiometries(
     known_value="cyclable lithium capacity",
     options=None,
     tol=1e-6,
+    inputs=None,
 ):
     """
     Calculate initial stoichiometries to start off the simulation at a particular
@@ -866,7 +871,7 @@ def get_initial_stoichiometries(
         The initial stoichiometries that give the desired initial state of charge
     """
     esoh_solver = ElectrodeSOHSolver(parameter_values, param, known_value, options)
-    return esoh_solver.get_initial_stoichiometries(initial_value, tol)
+    return esoh_solver.get_initial_stoichiometries(initial_value, tol, inputs=inputs)
 
 
 def get_min_max_stoichiometries(

pybamm/models/full_battery_models/lithium_ion/electrode_soh_half_cell.py

@@ -28,7 +28,7 @@ def __init__(self, name="ElectrodeSOH model"):
 
         x_100 = pybamm.Variable("x_100", bounds=(0, 1))
         x_0 = pybamm.Variable("x_0", bounds=(0, 1))
-        Q_w = pybamm.InputParameter("Q_w")
+        Q_w = pybamm.Parameter("Total lithium in the working electrode [mol]")
         T_ref = param.T_ref
         U_w = param.p.prim.U
         Q = Q_w * (x_100 - x_0)
@@ -64,6 +64,7 @@ def get_initial_stoichiometry_half_cell(
     param=None,
     known_value="cyclable lithium capacity",
     options=None,
+    inputs=None,
 ):
     """
     Calculate initial stoichiometry to start off the simulation at a particular
@@ -86,7 +87,7 @@ def get_initial_stoichiometry_half_cell(
         The initial stoichiometry that give the desired initial state of charge
     """
     param = pybamm.LithiumIonParameters(options)
-    x_0, x_100 = get_min_max_stoichiometries(parameter_values)
+    x_0, x_100 = get_min_max_stoichiometries(parameter_values, inputs=inputs)
 
     if isinstance(initial_value, str) and initial_value.endswith("V"):
         V_init = float(initial_value[:-1])
@@ -129,7 +130,7 @@ def get_initial_stoichiometry_half_cell(
     return x
 
 
-def get_min_max_stoichiometries(parameter_values, options=None):
+def get_min_max_stoichiometries(parameter_values, options=None, inputs=None):
     """
     Get the minimum and maximum stoichiometries from the parameter values
 
@@ -142,12 +143,20 @@ def get_min_max_stoichiometries(parameter_values, options=None):
         :class:`pybamm.BatteryModelOptions`.
         If None, the default is used: {"working electrode": "positive"}
     """
+    inputs = inputs or {}
     if options is None:
         options = {"working electrode": "positive"}
     esoh_model = pybamm.lithium_ion.ElectrodeSOHHalfCell("ElectrodeSOH")
     param = pybamm.LithiumIonParameters(options)
+    # Use Q_w as a symbol rather than a value
+    Q_w = param.p.Q_init
+    # Q_w = parameter_values.evaluate(param.p.Q_init)
+    # Add Q_w to parameter values
+    parameter_values.update(
+        {"Total lithium in the working electrode [mol]": Q_w},
+        check_already_exists=False,
+    )
     esoh_sim = pybamm.Simulation(esoh_model, parameter_values=parameter_values)
-    Q_w = parameter_values.evaluate(param.p.Q_init)
-    esoh_sol = esoh_sim.solve([0], inputs={"Q_w": Q_w})
+    esoh_sol = esoh_sim.solve([0], inputs=inputs)
     x_0, x_100 = esoh_sol["x_0"].data[0], esoh_sol["x_100"].data[0]
     return x_0, x_100

pybamm/parameters/parameter_values.py

@@ -292,14 +292,20 @@ def set_initial_stoichiometry_half_cell(
         known_value="cyclable lithium capacity",
         inplace=True,
         options=None,
+        inputs=None,
     ):
         """
         Set the initial stoichiometry of the working electrode, based on the initial
         SOC or voltage
         """
         param = param or pybamm.LithiumIonParameters(options)
         x = pybamm.lithium_ion.get_initial_stoichiometry_half_cell(
-            initial_value, self, param=param, known_value=known_value, options=options
+            initial_value,
+            self,
+            param=param,
+            known_value=known_value,
+            options=options,
+            inputs=inputs,
         )
         if inplace:
             parameter_values = self
@@ -324,6 +330,7 @@ def set_initial_stoichiometries(
         known_value="cyclable lithium capacity",
         inplace=True,
         options=None,
+        inputs=None,
         tol=1e-6,
     ):
         """
@@ -338,6 +345,7 @@ def set_initial_stoichiometries(
             known_value=known_value,
             options=options,
             tol=tol,
+            inputs=inputs,
         )
         if inplace:
             parameter_values = self
@@ -769,7 +777,7 @@ def _process_symbol(self, symbol):
             # Backup option: return the object
             return symbol
 
-    def evaluate(self, symbol):
+    def evaluate(self, symbol, inputs=None):
         """
         Process and evaluate a symbol.
 
@@ -787,7 +795,15 @@ def evaluate(self, symbol):
         if processed_symbol.is_constant():
             return processed_symbol.evaluate()
         else:
-            raise ValueError("symbol must evaluate to a constant scalar or array")
+            # In the case that the only issue is an input parameter contained in inputs,
+            # go ahead and try and evaluate it with the inputs. If it doesn't work, raise
+            # the value error.
+            try:
+                return processed_symbol.evaluate(inputs=inputs)
+            except Exception as exc:
+                raise ValueError(
+                    "symbol must evaluate to a constant scalar or array"
+                ) from exc
 
     def _ipython_key_completions_(self):
         return list(self._dict_items.keys())