diff --git a/changelog.d/fix-empty-group-reductions.fixed.md b/changelog.d/fix-empty-group-reductions.fixed.md
new file mode 100644
index 000000000..be0927fb2
--- /dev/null
+++ b/changelog.d/fix-empty-group-reductions.fixed.md
@@ -0,0 +1 @@
+Group `all`, `max` and `min`, `members_position` and `Population.get_rank` now work when a simulation has no persons, giving the neutral element (True, -inf, inf) for each group instead of raising "zero-size array to reduction"; memberships given as an empty list, and default simulations of no one, are kept as integers so ranks and projections over them work.
diff --git a/policyengine_core/populations/group_population.py b/policyengine_core/populations/group_population.py
index 70ff98125..8d42a811f 100644
--- a/policyengine_core/populations/group_population.py
+++ b/policyengine_core/populations/group_population.py
@@ -67,7 +67,11 @@ def clone(
     def members_position(self) -> ArrayLike:
         if self._members_position is None and self.members_entity_id is not None:
             # We could use self.count and self.members.count , but with the current initilization, we are not sure count will be set before members_position is called
-            nb_entities = numpy.max(self.members_entity_id) + 1
+            nb_entities = (
+                numpy.max(self.members_entity_id) + 1
+                if len(self.members_entity_id)
+                else 0
+            )
             nb_persons = len(self.members_entity_id)
             self._members_position = numpy.empty_like(self.members_entity_id)
             counter_by_entity = numpy.zeros(nb_entities)
@@ -88,6 +92,10 @@ def members_entity_id(self) -> ArrayLike:
 
     @members_entity_id.setter
     def members_entity_id(self, members_entity_id: ArrayLike) -> None:
+        if members_entity_id is not None and len(members_entity_id) == 0:
+            # An empty list or ``numpy.array([])`` is float, and NumPy indexes
+            # only with integers (or booleans): keep no members as integers.
+            members_entity_id = numpy.zeros(0, dtype=numpy.intp)
         self._members_entity_id = members_entity_id
 
     @property
@@ -191,7 +199,10 @@ def reduce(
 
         # We loop over the positions in the entity
         # Looping over the entities is tempting, but potentielly slow if there are a lot of entities
-        biggest_entity_size = numpy.max(position_in_entity) + 1
+        # With no persons there is no position, and every entity keeps the neutral element.
+        biggest_entity_size = (
+            numpy.max(position_in_entity) + 1 if len(position_in_entity) else 0
+        )
 
         for p in range(biggest_entity_size):
             values = self.value_nth_person(p, filtered_array, default=neutral_element)
diff --git a/policyengine_core/populations/population.py b/policyengine_core/populations/population.py
index d2e8df1da..fe0b37a03 100644
--- a/policyengine_core/populations/population.py
+++ b/policyengine_core/populations/population.py
@@ -233,17 +233,22 @@ def get_rank(
         )
 
         positions = entity.members_position
-        biggest_entity_size = numpy.max(positions) + 1
+        biggest_entity_size = numpy.max(positions) + 1 if len(positions) else 0
         filtered_criteria = numpy.where(condition, criteria, numpy.inf)
         ids = entity.members_entity_id
 
         # Matrix: the value in line i and column j is the value of criteria for the jth person of the ith entity
-        matrix = numpy.asarray(
-            [
-                entity.value_nth_person(k, filtered_criteria, default=numpy.inf)
-                for k in range(biggest_entity_size)
-            ]
-        ).transpose()
+        # (``reshape`` keeps it two-dimensional when there are no persons)
+        matrix = (
+            numpy.asarray(
+                [
+                    entity.value_nth_person(k, filtered_criteria, default=numpy.inf)
+                    for k in range(biggest_entity_size)
+                ]
+            )
+            .reshape(biggest_entity_size, entity.count)
+            .transpose()
+        )
 
         # We double-argsort all lines of the matrix.
         # Double-argsorting gets the rank of each value once sorted
diff --git a/policyengine_core/simulations/simulation_builder.py b/policyengine_core/simulations/simulation_builder.py
index e3dee0cc5..5f3840501 100644
--- a/policyengine_core/simulations/simulation_builder.py
+++ b/policyengine_core/simulations/simulation_builder.py
@@ -245,7 +245,7 @@ def build_default_simulation(
             )
         for population in simulation.populations.values():
             population.count = count
-            population.ids = np.array(range(count))
+            population.ids = np.arange(count)  # integers even when count is 0
             if not population.entity.is_person:
                 population.members_entity_id = (
                     population.ids
diff --git a/tests/core/test_group_population_empty.py b/tests/core/test_group_population_empty.py
new file mode 100644
index 000000000..60a6f005f
--- /dev/null
+++ b/tests/core/test_group_population_empty.py
@@ -0,0 +1,83 @@
+"""Group reductions and ranks when a simulation has no persons.
+
+``GroupPopulation.reduce`` (behind ``all``, ``max`` and ``min``) returns its
+neutral element for a group with no members, and ``sum`` and ``nb_persons``
+return 0. With no persons at all, ``reduce``, ``members_position`` and
+``Population.get_rank`` took ``numpy.max`` of an empty array and raised
+"zero-size array to reduction operation maximum which has no identity".
+Memberships given as an empty list, or built for no one, were also float,
+which NumPy cannot index with.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from policyengine_core.entities import build_entity
+from policyengine_core.populations import GroupPopulation, Population
+from policyengine_core.simulations import SimulationBuilder
+from policyengine_core.taxbenefitsystems import TaxBenefitSystem
+
+PERSON = build_entity("person", "persons", "Person", is_person=True)
+HOUSEHOLD = build_entity(
+    "household",
+    "households",
+    "Household",
+    roles=[{"key": "member", "plural": "members"}],
+)
+
+
+def populations(members_entity_id, count):
+    persons = Population(PERSON)
+    persons.count = len(members_entity_id)
+    households = GroupPopulation(HOUSEHOLD, persons)
+    households.count = count
+    households.members_entity_id = np.asarray(members_entity_id, dtype=int)
+    return persons, households
+
+
+def test_reductions_over_no_persons_give_neutral_elements():
+    _, households = populations([], 1)
+    no_values = np.array([], dtype=float)
+
+    assert households.all(no_values > 0).tolist() == [True]
+    assert households.any(no_values > 0).tolist() == [False]
+    assert households.max(no_values).tolist() == [-np.inf]
+    assert households.min(no_values).tolist() == [np.inf]
+    assert households.sum(no_values).tolist() == [0]
+    assert households.nb_persons().tolist() == [0]
+    assert households.members_position.tolist() == []
+
+
+def test_rank_over_no_persons_is_empty():
+    persons, households = populations([], 2)
+
+    assert persons.get_rank(households, np.array([], dtype=float)).tolist() == []
+
+
+def test_no_members_given_as_a_plain_empty_list():
+    # ``[]`` and ``numpy.array([])`` are float; NumPy indexes only with
+    # integers, so ranks and projections raised IndexError.
+    persons, households = populations([], 1)
+    households.members_entity_id = []
+
+    assert households.members_position.dtype.kind == "i"
+    assert persons.get_rank(households, np.array([], dtype=float)).tolist() == []
+    assert households.project(np.array([5.0])).tolist() == []
+    assert households.max(np.array([], dtype=float)).tolist() == [-np.inf]
+
+
+def test_default_simulation_of_no_one():
+    simulation = SimulationBuilder().build_default_simulation(
+        TaxBenefitSystem([PERSON, HOUSEHOLD]), count=0
+    )
+    persons = simulation.persons
+    households = simulation.populations["household"]
+    no_values = np.array([], dtype=float)
+
+    assert households.ids.dtype.kind == "i"
+    assert persons.get_rank(households, no_values).tolist() == []
+    assert households.value_nth_person(0, no_values).tolist() == []
+    assert households.project(no_values).tolist() == []
+    assert persons.household.max(no_values).tolist() == []
+    assert households.all(no_values > 0).tolist() == []
diff --git a/tests/core/test_group_population_reductions_property.py b/tests/core/test_group_population_reductions_property.py
new file mode 100644
index 000000000..2732cfc03
--- /dev/null
+++ b/tests/core/test_group_population_reductions_property.py
@@ -0,0 +1,80 @@
+"""Property: group reductions and ranks match a plain-Python oracle.
+
+For any number of groups (1 to 6), any assignment of 0 to 30 persons to them
+(so groups may be empty, and there may be no persons at all), given as an
+integer array or a plain list, and any values:
+
+* ``sum``, ``nb_persons``, ``any``, ``all``, ``max`` and ``min`` equal the
+  oracle, with 0, False, True, -inf and inf for groups with no members;
+* ``Population.get_rank`` gives each person the position of their value among
+  their group's values, ties going to the person listed first.
+
+Examples, including the empty-population witness, are in
+``test_group_population_empty.py``.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+# The smoke job installs Core without the dev extra but collects every module.
+pytest.importorskip("hypothesis")
+
+from hypothesis import given, settings  # noqa: E402
+from hypothesis import strategies as st  # noqa: E402
+
+from tests.core.test_group_population_empty import populations  # noqa: E402
+
+
+@st.composite
+def _groups(draw):
+    count = draw(st.integers(1, 6))
+    members_entity_id = draw(st.lists(st.integers(0, count - 1), max_size=30))
+    values = draw(
+        st.lists(
+            st.integers(-5, 5),
+            min_size=len(members_entity_id),
+            max_size=len(members_entity_id),
+        )
+    )
+    as_list = draw(st.booleans())
+    return count, members_entity_id, values, as_list
+
+
+@settings(max_examples=300, deadline=None)
+@given(_groups())
+def test_reductions_and_ranks_match_oracle(case):
+    count, members_entity_id, values, as_list = case
+    persons, households = populations(members_entity_id, count)
+    if as_list:
+        # Memberships may come as a plain list; an empty one is float.
+        households.members_entity_id = list(members_entity_id)
+    array = np.array(values, dtype=float)
+
+    members = {group: [] for group in range(count)}
+    for person, group in enumerate(members_entity_id):
+        members[group].append(person)
+    group_values = [[values[person] for person in members[g]] for g in range(count)]
+
+    assert households.sum(array).tolist() == [sum(v) for v in group_values]
+    assert households.nb_persons().tolist() == [len(v) for v in group_values]
+    assert households.any(array > 0).tolist() == [
+        any(x > 0 for x in v) for v in group_values
+    ]
+    assert households.all(array > 0).tolist() == [
+        all(x > 0 for x in v) for v in group_values
+    ]
+    assert households.max(array).tolist() == [
+        max(v, default=-np.inf) for v in group_values
+    ]
+    assert households.min(array).tolist() == [
+        min(v, default=np.inf) for v in group_values
+    ]
+
+    expected_rank = [0] * len(values)
+    for group in range(count):
+        ordered = sorted(members[group], key=lambda person: values[person])
+        for rank, person in enumerate(ordered):
+            expected_rank[person] = rank
+    assert persons.get_rank(households, array).tolist() == expected_rank