The detection of gravitational wave events has stimulated theoretical modeling of the formation and evolution of double compact objects (DCOs). However, even for the most studied isolated binary evolution channel, there exist large uncertainties in the input parameters and treatments of the binary evolution process. So far, double neutron stars (DNSs) are the only DCOs for which direct observations are available through traditional electromagnetic astronomy. In this work, we adopt a population synthesis method to investigate the formation and evolution of Galactic DNSs. We construct 324 models for the formation of Galactic DNSs, taking into account various possible combinations of critical input parameters and processes such as mass transfer efficiency, supernova type, common envelope efficiency, neutron star kick velocity, and pulsar selection effect. We employ Bayesian analysis to evaluate the adopted models by comparing with observations. We also compare the expected DNS merger rate in the galaxy with that inferred from the known Galactic population of pulsar-neutron star systems. Based on these analyses we derive the favorable range of the aforementioned key parameters.