I haven't seen anyone claim that BMS symmetries relate two physically inequivalent systems. It's important to realize that the "system" isn't just the part of the spacetime on one slice of future null infinity $\mathscr{I}^+$, or in one region of $\mathscr{I}^+$ — it's all of $\mathscr{I}^+$. The BMS group is originally defined as the group of all asymptotic gauge transformations that are symmetries of the asymptotic metric. In particular, any BMS transformation really is just a gauge transformation, which cannot affect the physical meaning of the quantities represented in terms of coordinates either before or after the BMS transformation.
But this is subtly different from what Hawking, Strominger, etc., are claiming. Instead, the claim really is about parts of the same physical system at two different times. The usual idea is that you can look at a very early time slice of $\mathscr{I}^+$, and another very late time slice — such that radiation was given off between the two times. [It's usually also stipulated that things are quiet (no radiation) before and after those two times, though it's not entirely clear to me how much that actually matters.] But the spacetime preserves some "memory" of the gravitational waves that were emitted. The way that the spacetime approaches future null infinity is actually different at early times and at late times.
Now, it's too much to explain all this here, so I'll just assume that you are familiar with the basic mathematics of gravitational-wave (GW) theory, and I'll gloss over some details. You'll also need to understand the BMS group — specifically the supertranslations. I think I did a decent job explaining them here, or in references therein. The key fact is that the retarded-time coordinate $u$ on $\mathscr{I}^+$ transforms under a supertranslation as
\begin{equation}
u \to u + \alpha(\theta, \phi),
\end{equation}
where $\alpha$ is an arbitrary (real-valued, and hopefully smooth) function of the spherical coordinates.
Basically, the GW information is carried in the metric perturbation $h$, which falls off as $1/r$ asymptotically. Obviously, the perturbation is precisely zero on $\mathscr{I}^+$, but we just talk about the leading-order term. You can extend the BMS transformation (which, strictly speaking, only applies on $\mathscr{I}^+$) to a transformation over an infinitesimal neighborhood of $\mathscr{I}^+$, and see how it affects the metric perturbation. It turns out that the metric perturbation transforms as
\begin{equation}
h \to h - \bar{\eth}^2 \alpha.
\end{equation}
(This is just the transformation under supertranslation, without any rotation or boost. Also, I've suppressed the complicated dependence on location.) Here, $\eth$ is an angular-derivative operator introduced by Newman and Penrose, and among other things, it removes any contribution from any time and space translations present in the supertranslation. Any BMS transformation, including the supertranslation $\alpha$, is constant in time. So this is just saying that you get to add a constant-in-time angle-dependent quantity to the metric perturbation. (Not an arbitrary constant-in-time angle-dependent quantity; remember that $\alpha$ must be real-valued.)
The claim, then, is that you can devise a supertranslation that lets you cancel out any GW memory effect on a given slice of $\mathscr{I}^+$. Now it would seem that this means that the memory effect is pure gauge and should be ignored. But that's not quite true. Instead, because the BMS transformation is constant in time, you can only cancel out the memory on one time slice — but if the memory changes ever again, you can't still cancel it. Or, to put it another way, you might say that while the memory on one time slice could be considered a gauge effect, the difference between memories on two slices is a physical effect.