diff --git a/repyability/non_repairable.py b/repyability/non_repairable.py
index 358ca4d..16b5588 100644
--- a/repyability/non_repairable.py
+++ b/repyability/non_repairable.py
@@ -2,9 +2,7 @@
 from scipy.integrate import quad
 from scipy.interpolate import interp1d
 from scipy.optimize import minimize
-from surpyval.nonparametric import NonParametric
-from surpyval.parametric import Parametric
-from surpyval.parametric.exact_event_time import ExactEventTime
+from surpyval import ExactEventTime, NonParametric, Parametric
 
 from repyability.rbd.standby_node import StandbyModel
 
@@ -63,6 +61,11 @@ def avg_replacement_time(self, t):
         return out
 
     def _cost_rate(self, t):
+        # surpyval >= 0.11 returns 1-element arrays from sf/mean/qf, and t can
+        # arrive as a 1-element array from scipy.optimize.minimize. Coerce to a
+        # plain float so the quadrature bound in avg_replacement_time() is a
+        # scalar (scipy.integrate.quad rejects array bounds).
+        t = float(np.asarray(t).item())
         planned_costs = self.reliability_function(t) * self.cp
         unplanned_costs = (1 - self.reliability_function(t)) * self.cu
         avg_repl_time = self.avg_replacement_time(t)
diff --git a/repyability/rbd/helper_classes.py b/repyability/rbd/helper_classes.py
index 5741c1b..d5cd23b 100644
--- a/repyability/rbd/helper_classes.py
+++ b/repyability/rbd/helper_classes.py
@@ -10,6 +10,11 @@ def sf(cls, x):
     def ff(cls, x):
         return np.zeros_like(x).astype(float)
 
+    @classmethod
+    def cs(cls, x, X):
+        # Always working -> conditional survival is always 1.
+        return np.ones_like(np.atleast_1d(x)).astype(float)
+
     @classmethod
     def random(cls, size):
         return np.ones(size) * np.inf
@@ -24,6 +29,11 @@ def sf(cls, x):
     def ff(cls, x):
         return np.ones_like(x).astype(float)
 
+    @classmethod
+    def cs(cls, x, X):
+        # Never working -> conditional survival is always 0.
+        return np.zeros_like(np.atleast_1d(x)).astype(float)
+
     @classmethod
     def random(self, size):
         return np.zeros(size)
diff --git a/repyability/rbd/non_repairable_rbd.py b/repyability/rbd/non_repairable_rbd.py
index 439d2f7..bb4491e 100644
--- a/repyability/rbd/non_repairable_rbd.py
+++ b/repyability/rbd/non_repairable_rbd.py
@@ -10,6 +10,8 @@
 from numpy.typing import ArrayLike
 from surpyval import NonParametric
 
+from repyability.utils.wrappers import conditional_survival
+
 from .helper_classes import PerfectReliability, PerfectUnreliability
 from .rbd import RBD
 from .repeated_node import RepeatedNode
@@ -327,6 +329,28 @@ def reliability(self, x: Optional[ArrayLike] = None, *args, **kwargs):
         """
         return self.sf(x, *args, **kwargs)
 
+    def cs(self, x: ArrayLike, X: ArrayLike, *args, **kwargs) -> np.ndarray:
+        """Returns the conditional survival of the system.
+
+        That is, the probability the system survives a *further* ``x`` given it
+        has already survived to ``X``: ``R(x | X) = sf(X + x) / sf(X)``.
+
+        Parameters
+        ----------
+        x : ArrayLike
+            The further duration/s at which conditional survival is evaluated.
+        X : ArrayLike
+            The age/s the system is known to have survived to.
+        *args, **kwargs :
+            Any sf() arguments (e.g. working_nodes, broken_nodes, method).
+
+        Returns
+        -------
+        np.ndarray
+            The conditional survival probability/ies.
+        """
+        return conditional_survival(self, x, X, *args, **kwargs)
+
     @check_x
     def sf_by_node(
         self, x: Optional[ArrayLike] = None, *args, **kwargs
@@ -340,11 +364,10 @@ def sf_by_node(
             node_sf[node_name] = node.sf(x)
         return node_sf
 
+    @check_x
     def ff_by_node(
         self, x: Optional[ArrayLike] = None, *args, **kwargs
     ) -> Dict[Any, np.ndarray]:
-
-        # The return dict
         node_ff: Dict[Any, np.ndarray] = {}
 
         # Cache the component reliabilities for efficiency
@@ -437,7 +460,11 @@ def random(self, size):
         for i in range(size):
             event_queue: PriorityQueue = PriorityQueue()
             for node in self.G.nodes:
-                time = self.reliabilities[node].random(1)
+                # .random(1) returns a 1-element array; take the scalar so the
+                # event time orders the PriorityQueue and assigns into ``out``
+                # (NumPy >= 2 rejects assigning a 1-element array to a scalar).
+                draw = np.asarray(self.reliabilities[node].random(1))
+                time = float(draw.reshape(-1)[0])
                 event_queue.put(NodeFailure(time, node))
 
             working_nodes = {k: True for k in self.G.nodes}
diff --git a/repyability/rbd/repeated_node.py b/repyability/rbd/repeated_node.py
index 14ff116..74a8c8c 100644
--- a/repyability/rbd/repeated_node.py
+++ b/repyability/rbd/repeated_node.py
@@ -44,3 +44,9 @@ def ff(self, x):
         else:
             ff = self.model.ff(x) ** self.repeats
         return ff
+
+    def cs(self, x, X):
+        """Conditional survival ``R(x | X) = sf(X + x) / sf(X)``."""
+        from repyability.utils.wrappers import conditional_survival
+
+        return conditional_survival(self, x, X)
diff --git a/repyability/rbd/repeated_standby_node.py b/repyability/rbd/repeated_standby_node.py
index 4b38825..ce833ac 100644
--- a/repyability/rbd/repeated_standby_node.py
+++ b/repyability/rbd/repeated_standby_node.py
@@ -60,3 +60,9 @@ def sf(self, *args, **kwargs):
 
     def ff(self, *args, **kwargs):
         return self._sf_model.ff(*args, **kwargs)
+
+    def cs(self, x, X):
+        """Conditional survival ``R(x | X) = sf(X + x) / sf(X)``."""
+        from repyability.utils.wrappers import conditional_survival
+
+        return conditional_survival(self, x, X)
diff --git a/repyability/rbd/standby_node.py b/repyability/rbd/standby_node.py
index f96eb18..613ba1f 100644
--- a/repyability/rbd/standby_node.py
+++ b/repyability/rbd/standby_node.py
@@ -179,3 +179,9 @@ def ff(self, *args, **kwargs):
         if self._sf_model is not None:
             return self._sf_model.ff(*args, **kwargs)
         return self.model.ff(*args, **kwargs)
+
+    def cs(self, x, X):
+        """Conditional survival ``R(x | X) = sf(X + x) / sf(X)``."""
+        from repyability.utils.wrappers import conditional_survival
+
+        return conditional_survival(self, x, X)
diff --git a/repyability/tests/test_rbd_misc.py b/repyability/tests/test_rbd_misc.py
index 72d92b7..7f456e2 100644
--- a/repyability/tests/test_rbd_misc.py
+++ b/repyability/tests/test_rbd_misc.py
@@ -3,12 +3,28 @@
 
 Uses pytest fixtures located in conftest.py in the tests/ directory.
 """
+import numpy as np
 import pytest
-from surpyval import FixedEventProbability
+from surpyval import Exponential, FixedEventProbability
 
 from repyability.rbd.non_repairable_rbd import NonRepairableRBD
 
 
+def test_rbd_mean_time_to_failure_series():
+    # Two exponential components in series fail at the first failure, so the
+    # system MTTF is 1 / (lambda_1 + lambda_2). Exercises the Monte-Carlo
+    # random()/mean() path.
+    rbd = NonRepairableRBD(
+        [("input", "a"), ("a", "b"), ("b", "output")],
+        {
+            "a": Exponential.from_params([0.01]),
+            "b": Exponential.from_params([0.02]),
+        },
+    )
+    np.random.seed(0)
+    assert rbd.mean(mc_samples=20_000) == pytest.approx(1 / 0.03, rel=5e-2)
+
+
 # Check components are correct lengths
 def test_rbd_components(rbd1: NonRepairableRBD, rbd2: NonRepairableRBD):
     assert len(rbd1.nodes) == 3
diff --git a/repyability/tests/test_rbd_sf.py b/repyability/tests/test_rbd_sf.py
index 76b35db..bcfb77a 100644
--- a/repyability/tests/test_rbd_sf.py
+++ b/repyability/tests/test_rbd_sf.py
@@ -15,7 +15,9 @@ def test_rbd_sf_ff_by_node(rbd_series: NonRepairableRBD):
     node_sf = rbd_series.sf_by_node(t)
     node_ff = rbd_series.ff_by_node(t)
     for k in node_sf.keys():
-        assert node_ff[k] == 1 - node_sf[k]
+        # surpyval >= 0.11 computes ff(t) directly (more precisely than
+        # 1 - sf(t)), so the two can differ by a floating-point ULP.
+        assert node_ff[k] == pytest.approx(1 - node_sf[k])
 
 
 def test_rel_unrel(rbd_series: NonRepairableRBD):
@@ -144,3 +146,23 @@ def test_rbd_sf_NonRepairableRBD_as_node(rbd1: NonRepairableRBD, method: str):
         * rbd.reliabilities["NonRepairableRBD"].ff(t)
         * rbd.reliabilities[4].ff(t)
     ) == rbd.ff(t, method=method)
+
+
+def test_rbd_conditional_survival():
+    """cs(x, X) returns sf(X + x) / sf(X), clipped to [0, 1]."""
+    import numpy as np
+    from surpyval import Exponential, Weibull
+
+    rbd = NonRepairableRBD(
+        [("input", "c1"), ("c1", "c2"), ("c2", "output")],
+        {
+            "c1": Weibull.from_params([100, 2]),
+            "c2": Exponential.from_params([0.01]),
+        },
+    )
+    x = np.array([10.0, 20.0, 50.0])
+    X = 30.0
+    expected = rbd.sf(x + X) / rbd.sf(X)
+    assert np.allclose(rbd.cs(x, X), expected)
+    # Surviving zero further time is certain.
+    assert rbd.cs(0.0, X)[0] == 1.0
diff --git a/repyability/utils/wrappers.py b/repyability/utils/wrappers.py
index cdf4c54..21024c2 100644
--- a/repyability/utils/wrappers.py
+++ b/repyability/utils/wrappers.py
@@ -1,3 +1,6 @@
+import numpy as np
+
+
 def check_probability(func):
     """checks probability is between 0 and 1"""
 
@@ -10,3 +13,38 @@ def wrap(obj, target: float, *args, **kwargs):
             return func(obj, target, *args, **kwargs)
 
     return wrap
+
+
+def conditional_survival(model, x, X, *args, **kwargs):
+    """Conditional survival from any model exposing ``sf``.
+
+    Returns the probability of surviving a *further* ``x`` given the item has
+    already survived to ``X``:
+
+    .. math::
+        R(x \\mid X) = \\frac{R(X + x)}{R(X)}
+
+    Parameters
+    ----------
+    model : object
+        Anything with an ``sf(x, ...)`` method (a distribution, a standby
+        arrangement, an RBD, ...).
+    x : array_like or scalar
+        The further duration(s) at which conditional survival is evaluated.
+    X : array_like or scalar
+        The age(s) the item is known to have survived to.
+
+    Returns
+    -------
+    numpy.ndarray
+        The conditional survival probability, clipped to ``[0, 1]``; where the
+        item has all but surely failed by ``X`` (``R(X) ≈ 0``) it is ``0``.
+    """
+    x = np.atleast_1d(np.asarray(x, dtype=float))
+    X = np.atleast_1d(np.asarray(X, dtype=float))
+    denom = np.asarray(model.sf(X, *args, **kwargs), dtype=float)
+    numer = np.asarray(model.sf(x + X, *args, **kwargs), dtype=float)
+    with np.errstate(divide="ignore", invalid="ignore"):
+        out = numer / denom
+    out = np.where(np.isfinite(out), out, 0.0)
+    return np.clip(out, 0.0, 1.0)
diff --git a/requirements.txt b/requirements.txt
index 6ef5fd1..6209e49 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,6 +1,6 @@
 # Production requirements
 # For development, use requirements_dev.txt instead
-networkx==3.0
-numpy==1.23.5
-surpyval==0.10.10
-tqdm==4.64.1
+networkx>=3.0
+numpy>=2.0,<3
+surpyval==0.11.1
+tqdm>=4.64