  We propose and construct a two-parameter perturbative expansion around a
Friedmann-Lema\^{i}tre-Robertson-Walker geometry that can be used to model
high-order gravitational effects in the presence of non-linear structure. This
framework reduces to the weak-field and slow-motion post-Newtonian treatment of
gravity in the appropriate limits, but also includes the low-amplitude
large-scale fluctuations that are important for cosmological modelling. We
derive a set of field equations that can be applied to the late Universe, where
non-linear structure exists on supercluster scales, and perform a detailed
investigation of the associated gauge problem. This allows us to identify a
consistent set of perturbed quantities in both the gravitational and matter
sectors, and to construct a set of gauge-invariant quantities that correspond
to each of them. The field equations, written in terms of these quantities,
take on a relatively simple form, and allow the effects of small-scale
structure on the large-scale properties of the Universe to be clearly
identified. We find that inhomogeneous structures source the global expansion,
that there exist new field equations at new orders, and that there is vector
gravitational potential that is a hundred times larger than one might naively
expect from cosmological perturbation theory. Finally, we expect our formalism
to be of use for calculating relativistic effects in upcoming ultra-large-scale
surveys, as the form of the gravitational coupling between small and large
scales depends on the non-linearity of Einstein's equations, and occurs at what
is normally thought of as first order in cosmological perturbations.
