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 electromagnetic observations are available. In this work, we adopt a population synthesis method to investigate the formation and evolution of Galactic DNSs. We construct 216 models for the formation of Galactic DNSs, taking into account various possible combinations of critical input parameters and processes such as mass transfer (MT) efficiency, supernova type, common envelope (CE) 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 inferrd from the known Galactic population of Pulsar-NS systems. Based on these analyses we derive favorable range of the aforementioned key parameters.