Abstract Kilonovae are optical flashes produced in the aftermath of neutron star-neutron star mergers or neutron star-black hole mergers. In this work, we use the Millennium Simulation, combined with semi-analytic galaxy formation model GABE (Galaxy Assembly with Binary Evolution) to explore the cosmic event rate of kilonovae, and the properties of their host galaxies in a cosmological context. We find that model with supernova kick velocity of Vkick = 0 km s−1 fits the observation best, satisfying the preference for low kick velocity binary system in theoretical models. With Vkick = 0 km s−1, the cosmic event rate of NNMs and NBMs at z = 0 are 283 and 91 Gpc−3 yr−1, respectively, marginally consistent with the constraint from LVC GWTC-1. For Milky Way-mass galaxies, we predict the NNM rate is $25.7^{+59.6}_{-7.1}\, {\rm Myr}^{-1}$, which is also in good agreement with the observed properties of binary neutron stars in the Milky Way. Taking all the kilonovae into account in the history of Milky Way-mass galaxies, we find that the averaged r-process elements yield (A > 79) in a NNM/NBM event should be 0.01 M⊙ to be consistent with observation. We conclude that NGC 4993, the host galaxy of GW170817, is a typical host galaxy for NNMs. However, generally NNMs and NBMs tend to reside in young, blue, star-forming, late-type galaxies, with stellar mass and gaseous metallicity distribution peaking at M* = 1010.65 M⊙ and 12 + log (O/H) = 8.72 − 8.85, respectively. By exploring the connection between kilonovae and their host galaxies in the cosmological background, it is promising to constrain model better when we have more events in the forthcoming future.