Kilonovae are optical flashes produced in the aftermath of neutron star-neutron star mergers (NNMs) or neutron star-black hole mergers (NBMs). The multi-messager observation of the recent gravitational wave event GW170817 confirms that it originated from a NNM and triggered a kilonova. In this work, we use the Millennium Simulation, combined with a semianalytic galaxy formation model-GABE (Galaxy Assembly with Binary Evolution) which adopts binary stellar population synthesis models, 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 V kick = 0 km s −1 fits the observation best, in agreement with the exception of some formation channels of binary neutron star. This indicates that NNMs prefer to originate from binary systems with low kick velocities. With V kick = 0 km s −1 , the cosmic event rate of NNMs and NBMs at z = 0 are 283 Gpc −3 yr −1 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 Myr −1 , which is also in good agreement with the observed properties of binary neutron stars in the Milky Way. Taking all the NNMs into account in the history of Milky Way-mass galaxies, we find that the averaged r-process elements yield with A > 79 in a NNM and 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 * = 10 10.65 M ⊙ and 12 + log (O/H) = 8.72 − 8.85, respectively. By studying kilonovae host galaxies in the cosmological background, it is promising to constrain model details better when we have more events in the forthcoming future.