We predict the rate for exclusive double-charmonium production in electron-positron annihilation $e^+ e^- \to J/\psi+\eta_c$ using perturbative quantum chromodynamics and the NRQCD framework for hard, heavy-quarkonium exclusive processes. The cross sections measured at the $B$-factories Belle and Babar at $\sqrt{s}=10.6$ GeV disagree with the pQCD leading-order predictions by an order of magnitude. The predictions at next-to-leading order are, however, very sensitive to the choice of the renormalization scale, resulting in an apparent discrepancy between the theoretical prediction and the data. We show that this discrepancy can in fact be eliminated by applying the Principle of Maximum Conformality (PMC) to set the renormalization scale. ... By carefully applying the PMC to different topologies of the annihilation process, one achieves precise pQCD predictions, together with improved perturbative convergence. We also observe that the single-photon-fragmentation QED correction is important, an effect which increases the total cross-section by about $10\%$. The scale-fixed, scheme-independent cross-section predicted by the PMC is $\sigma_{\rm tot}|_{\rm PMC}=20.35 ^{+3.5}_{-3.8}$ fb, where the uncertainties come from the squared average of the errors due to the value of the charm mass and the uncertainty from the quarkonium wavefunctions at the origin. We find that the typical momentum flow of the process is $2.30$ GeV, which explains the guessed choice of $2-3$ GeV using conventional scale-setting. The scale-fixed $e^+ e^- \to J/\psi+\eta_c$ cross-section predicted by the PMC shows excellent agreement with the Belle and Babar measurements, emphasizing the importance of a rigorous renormalization scale-setting procedure. Show Partial Abstract