  Kilonovae are optical flashes produced in the aftermath of neutron
star-neutron star mergers (NNMs) or neutron star-black hole mergers (NBMs). In
this work, we use the Millennium Simulation, combined with a semi-analytic
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 0 km/s 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$_{\rm kick}$=0 km/s, 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$_{\odot}$ 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
$10^{10.65}$ M$_{\odot}$ and 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.
(abridged)
