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 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 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 } _ { \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 M _ { * } = 10 ^ { 10.65 } { M } _ { \odot } 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 .