  In the advent of large scale surveys, individually modelling
strong-gravitational lenses and their counterpart time-delays in order to
precisely measure $H_0$ will become computationally expensive, and highly
complex. A complimentary approach is to study the cumulative distribution
function (CDF) of time-delays where the global population of lenses is modelled
along with $H_0$. In this paper we use a suite of hydro-dynamical simulations
to estimate the CDF of time-delays from doubly-imaged quasars for a realistic
distribution of lenses. We find that the CDFs exhibit large amounts of
halo-halo variance, regulated by the density profile inner slope and the total
mass within $5$kpc. With the objective of fitting to data, we compress the CDFs
using Principal Component Analysis and fit a Gaussian Processes Regressor
consisting of three physical features: the redshift of the lens, $z_{\rm L}$;
the power law index of the halo, $\alpha$, and the mass within $5$kpc, plus
four cosmological features. Assuming a flat Universe, we fit our model to 27
doubly-imaged quasars finding $H_0=71^{+2}_{-3}$km/s/Mpc, $z_{\rm L} =
0.36_{-0.09}^{+0.2} $, $\alpha=-1.8_{-0.1}^{+0.1}$,
$\log(M(<5$kpc$)/M_\odot)=11.1_{-0.1}^{+0.1} $, $\Omega_{\rm M} =
0.3_{-0.04}^{+0.04} $ and $\Omega_{\rm \Lambda}=0.7_{-0.04}^{+0.04}$. We
compare our estimates of $z_{\rm L}$ and $\log(M(<5$kpc$)/M_\odot)$ to the data
and find that within the sensitivity of the data, they are not systematically
biased. We generate mock CDFs and find with that the Vera Rubin Observatory
(VRO) could measure $\sigma/H_0$ to $<3\%$, limited by the precision of the
model. If we are to exploit fully VRO, we require simulations that sample a
larger proportion of the lens population, with a variety of feedback models,
exploring all possible systematics.
