Utilizing first-order beam dynamics models is adequate for studying the beam properties during the conceptual design of a cyclotron-based proton therapy beamline. After finishing lattice design, particle-matter interaction simulations for passive elements (e.g., degrader, collimators, energy slit) are required. The cascade simulation is used for lattice updates in each iteration, which is complicated. In addition, when the models involve particle tracking and particle-matter interaction, their optimization process is time-consuming. Therefore, this study proposes a start-to-end modeling method using Monte Carlo Beam Delivery Simulation (BDSIM) software that considers more realistic factors, such as particle-matter interaction and the realistic vacuum chamber, to precisely evaluate working parameters, along with an efficient optimization method that utilizes multi-objective Bayesian optimization (MOBO) to improve transmission efficiency. Taking the Huazhong University of Science and Technology proton therapy facility (HUST-PTF) as an example, beam loss along the beamline is located, quantified, and subsequently reduced by tuning the quadrupole strengths based on MOBO. The results show that: (i) By considering the particle-matter interaction and the realistic vacuum chamber, the precision in the prediction of the beam properties is improved; (ii) After optimization, the transmission efficiency of the entire beamline is relatively increased by an average of 6.52 % under different energy settings, especially 11.39 % at 70 MeV.