Fix experiment initial guess regression

CHANGELOG.md

@@ -12,6 +12,8 @@
 
 ## Bug fixes
 
+- Fixed a regression causing inaccurate initial guesses at experiment step transitions. ([#5452](https://github.com/pybamm-team/PyBaMM/pull/5452))
+
 ## Breaking changes
 
 # [v26.3.1](https://github.com/pybamm-team/PyBaMM/tree/v26.3.1) - 2026-04-10

src/pybamm/models/base_model.py

@@ -1495,11 +1495,16 @@ def _resolve_from_var(target_var):
         equations = []
         for _target_slice, target_var, from_var, from_slice in entries:
             if from_var is None:
-                # Keep the target-model initial condition for unmapped variables
-                physical = (
-                    target_var.reference
-                    + target_var.scale * self.initial_conditions[target_var]
-                )
+                # If the variable is not a state variable of the previous model,
+                # query it by name. See #5449 for more details
+                try:
+                    physical = from_model.get_processed_variable(target_var.name)
+                except KeyError as e:
+                    raise pybamm.ModelError(
+                        f"Cannot build initial state mapper: variable "
+                        f"'{target_var.name}' is not a state in the previous "
+                        f"model and has no matching named variable to reuse."
+                    ) from e
             else:
                 physical = from_var.reference + from_var.scale * pybamm.StateVector(
                     from_slice

tests/integration/test_simulations.py

@@ -1,3 +1,4 @@
+import casadi
 import numpy as np
 import pytest
 
@@ -52,3 +53,60 @@ def test_pickle_roundtrip_exact(
                 np.testing.assert_array_equal(
                     sol_orig.yp, sol_loaded.yp, err_msg=f"{tag} sol.yp mismatch"
                 )
+
+
+class TestSimulationConsistentState:
+    def test_cv_initial_guess_uses_previous_current(self):
+        charge_current = 15.0
+        experiment = pybamm.Experiment(
+            [
+                (
+                    f"Charge at {charge_current}A until 4.2 V",
+                    "Hold at 4.2 V for 2 seconds",
+                )
+            ]
+        )
+        sim = pybamm.Simulation(
+            pybamm.lithium_ion.DFN(),
+            parameter_values=pybamm.ParameterValues("Chen2020"),
+            experiment=experiment,
+        )
+        sol = sim.solve(initial_soc=0.05)
+
+        # The mapper only runs when the two steps build distinct models.
+        cc_step, cv_step = sim.experiment.steps
+        cc_model = sim.steps_to_built_models[cc_step.basic_repr()]
+        cv_model = sim.steps_to_built_models[cv_step.basic_repr()]
+        assert cc_model is not cv_model
+
+        cc_solution = sol.cycles[0].steps[0]
+
+        # Replay BaseSolver.step's mapper invocation exactly: pull the
+        # compiled (func, jac_y, jac_p) tuple, build the CV step's input
+        # vector through the same code path, and call the func on the CC
+        # step's last raw state vector.
+        mapper_func, _, _ = sim._compiled_model_state_mappers[(cc_model, cv_model)]
+        cv_inputs_dict = sim._build_experiment_step_inputs(
+            {},
+            cv_step,
+            float(cc_solution.t[-1]),
+            None,
+            include_temperature=False,
+        )
+        cv_solver = sim._get_built_experiment_solver(cv_step)
+        model_inputs = cv_solver._set_up_model_inputs(cv_model, cv_inputs_dict)
+        p_vec = casadi.vertcat(*model_inputs.values())
+        y_from = cc_solution.last_state.all_ys[0]
+        seed = np.asarray(mapper_func(float(cc_solution.t[-1]), y_from, p_vec)).ravel()
+
+        # Decode the CV model's "Current variable [A]" entry from the seed.
+        current_var = next(
+            v for v in cv_model.y_slices if v.name == "Current variable [A]"
+        )
+        current_slice = cv_model.y_slices[current_var][0]
+        seed_current_scaled = float(seed[current_slice][0])
+        seed_current_phys = float(current_var.reference.evaluate()) + (
+            float(current_var.scale.evaluate()) * seed_current_scaled
+        )
+
+        assert seed_current_phys == pytest.approx(-charge_current, rel=1e-9)

tests/unit/test_models/test_base_model.py

@@ -1488,6 +1488,120 @@ def test_set_initial_conditions_symbol_with_value_attribute(self):
         disc_var2 = next(iter(model2_disc.initial_conditions.keys()))
         assert isinstance(model2_disc.initial_conditions[disc_var2], pybamm.Vector)
 
+    def test_build_initial_state_mapper_reuses_input_value(self):
+        # ``from_model``: only ``s`` is a state. ``current`` is exposed as a
+        # named variable whose value comes from an InputParameter.
+        s_from = pybamm.Variable("s")
+        from_model = pybamm.BaseModel()
+        from_model.rhs = {s_from: -s_from}
+        from_model.initial_conditions = {s_from: pybamm.Scalar(0.5)}
+        from_model.variables = {
+            "s": s_from,
+            "current": pybamm.InputParameter("I_in"),
+        }
+        pybamm.Discretisation().process_model(from_model)
+
+        # ``to_model``: ``current`` is now an algebraic state with a *wrong*
+        # initial-condition guess of 0 (must NOT be used by the mapper).
+        s_to = pybamm.Variable("s")
+        c_to = pybamm.Variable("current")
+        to_model = pybamm.BaseModel()
+        to_model.rhs = {s_to: -s_to}
+        to_model.algebraic = {c_to: c_to - pybamm.Scalar(1.0)}
+        to_model.initial_conditions = {
+            s_to: pybamm.Scalar(0.0),
+            c_to: pybamm.Scalar(0.0),
+        }
+        to_model.variables = {"s": s_to, "current": c_to}
+        pybamm.Discretisation().process_model(to_model)
+
+        mapper = to_model.build_initial_state_mapper(from_model)
+
+        # Evaluate the mapper with a known previous state vector and input.
+        y_from = np.array([[0.7]])
+        result = np.asarray(
+            mapper.evaluate(t=0, y=y_from, inputs={"I_in": 3.0})
+        ).ravel()
+
+        # Map (s, current) onto y_slices ordering of the target model.
+        s_slice = to_model.y_slices[s_to][0]
+        c_slice = to_model.y_slices[c_to][0]
+        assert result[s_slice][0] == pytest.approx(0.7)
+        # Critical: must equal the previous step's input, not the target IC (0).
+        assert result[c_slice][0] == pytest.approx(3.0)
+
+    def test_build_initial_state_mapper_handles_scale_and_reference(self):
+        from_scale_s = pybamm.Scalar(2.0)
+        from_ref_s = pybamm.Scalar(10.0)
+        s_from = pybamm.Variable("s", scale=from_scale_s, reference=from_ref_s)
+        from_model = pybamm.BaseModel()
+        from_model.rhs = {s_from: -s_from}
+        from_model.initial_conditions = {s_from: pybamm.Scalar(0.0)}
+        from_model.variables = {
+            "s": s_from,
+            "current": pybamm.InputParameter("I_in"),
+        }
+        pybamm.Discretisation().process_model(from_model)
+
+        to_scale_s = pybamm.Scalar(4.0)
+        to_ref_s = pybamm.Scalar(-5.0)
+        to_scale_c = pybamm.Scalar(0.5)
+        to_ref_c = pybamm.Scalar(100.0)
+        s_to = pybamm.Variable("s", scale=to_scale_s, reference=to_ref_s)
+        c_to = pybamm.Variable("current", scale=to_scale_c, reference=to_ref_c)
+        to_model = pybamm.BaseModel()
+        to_model.rhs = {s_to: -s_to}
+        to_model.algebraic = {c_to: c_to - pybamm.Scalar(1.0)}
+        to_model.initial_conditions = {
+            s_to: pybamm.Scalar(0.0),
+            c_to: pybamm.Scalar(0.0),
+        }
+        to_model.variables = {"s": s_to, "current": c_to}
+        pybamm.Discretisation().process_model(to_model)
+
+        mapper = to_model.build_initial_state_mapper(from_model)
+
+        # Choose physical values, encode via from_model scaling, decode in mapper.
+        phys_s = 17.0
+        y_s_from = (phys_s - 10.0) / 2.0  # 3.5
+        phys_current = 6.25
+        y_from = np.array([[y_s_from]])
+        result = np.asarray(
+            mapper.evaluate(t=0, y=y_from, inputs={"I_in": phys_current})
+        ).ravel()
+
+        s_slice = to_model.y_slices[s_to][0]
+        c_slice = to_model.y_slices[c_to][0]
+        # s: matched-state branch must round-trip through physical units.
+        # expected = (17 - (-5)) / 4 = 5.5
+        assert result[s_slice][0] == pytest.approx((phys_s - (-5.0)) / 4.0)
+        # current: input-parameter branch.
+        # expected = (6.25 - 100) / 0.5 = -187.5
+        assert result[c_slice][0] == pytest.approx((phys_current - 100.0) / 0.5)
+
+    def test_build_initial_state_mapper_missing_variable_raises(self):
+        s_from = pybamm.Variable("s")
+        from_model = pybamm.BaseModel()
+        from_model.rhs = {s_from: -s_from}
+        from_model.initial_conditions = {s_from: pybamm.Scalar(0.5)}
+        from_model.variables = {"s": s_from}
+        pybamm.Discretisation().process_model(from_model)
+
+        s_to = pybamm.Variable("s")
+        new_var = pybamm.Variable("brand_new")
+        to_model = pybamm.BaseModel()
+        to_model.rhs = {s_to: -s_to}
+        to_model.algebraic = {new_var: new_var - pybamm.Scalar(1.0)}
+        to_model.initial_conditions = {
+            s_to: pybamm.Scalar(0.0),
+            new_var: pybamm.Scalar(0.0),
+        }
+        to_model.variables = {"s": s_to, "brand_new": new_var}
+        pybamm.Discretisation().process_model(to_model)
+
+        with pytest.raises(pybamm.ModelError, match="brand_new"):
+            to_model.build_initial_state_mapper(from_model)
+
     def test_set_variables_error(self):
         var = pybamm.Variable("var")
         model = pybamm.BaseModel()