diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2eea9ed7c0..12c6ad8c84 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## Bug fixes
 
+- Fixed a bug where simulations using the CasADi-based solvers would fail randomly with the half-cell model ([#3494](https://github.com/pybamm-team/PyBaMM/pull/3494))
 - Fixed bug that made identical Experiment steps with different end times crash ([#3516](https://github.com/pybamm-team/PyBaMM/pull/3516))
 - Fixed bug in calculation of theoretical energy that made it very slow ([#3506](https://github.com/pybamm-team/PyBaMM/pull/3506))
 - The irreversible plating model now increments `f"{Domain} dead lithium concentration [mol.m-3]"`, not `f"{Domain} lithium plating concentration [mol.m-3]"` as it did previously. ([#3485](https://github.com/pybamm-team/PyBaMM/pull/3485))
diff --git a/pybamm/simulation.py b/pybamm/simulation.py
index f9aebb1c54..f743f4fc0f 100644
--- a/pybamm/simulation.py
+++ b/pybamm/simulation.py
@@ -4,6 +4,7 @@
 import pickle
 import pybamm
 import numpy as np
+import hashlib
 import warnings
 import sys
 from functools import lru_cache
@@ -133,6 +134,9 @@ def __init__(
         self._solution = None
         self.quick_plot = None
 
+        # Initialise instances of Simulation class with the same random seed
+        self._set_random_seed()
+
         # ignore runtime warnings in notebooks
         if is_notebook():  # pragma: no cover
             import warnings
@@ -156,6 +160,18 @@ def __setstate__(self, state):
         self.__dict__ = state
         self.get_esoh_solver = lru_cache()(self._get_esoh_solver)
 
+    # If the solver being used is CasadiSolver or its variants, set a fixed
+    # random seed during class initialization to the SHA-256 hash of the class
+    # name for purposes of reproducibility.
+    def _set_random_seed(self):
+        if isinstance(self._solver, pybamm.CasadiSolver) or isinstance(
+            self._solver, pybamm.CasadiAlgebraicSolver
+        ):
+            np.random.seed(
+                int(hashlib.sha256(self.__class__.__name__.encode()).hexdigest(), 16)
+                % (2**32)
+            )
+
     def set_up_and_parameterise_experiment(self):
         """
         Set up a simulation to run with an experiment. This creates a dictionary of