Vertical thermosyphon reboilers (VTRs) are widely used in distillation due to their high heat transfer efficiency and low fouling tendency. However, conventional design approaches and commercial tools often treat two-phase heat transfer, geometry, and pressure balance in isolation and depend on extensive manual tuning, which may yield low solution quality. To address this gap, a rigorous mixed integer nonlinear programming (MINLP) optimization framework is established for VTR design that couples thermodynamic, hydraulic, and geometric decisions simultaneously. A two-stage heat transfer modelling approach is proposed to capture the convective boiling transition by integrating sensible heating and boiling heat transfer mechanisms. Within this approach, a dynamic switching scheme selects heat transfer correlations according to the vapor fraction, which improves accuracy across different operating conditions. Moreover, a comprehensive pressure balance model is established that spans the external piping, the reboiler, and the column sump liquid level to ensure the resulting natural circulation is stable. Geometric and operating variables are optimized simultaneously to deliver coordinated design choices under duty and layout constraints. The optimization is implemented in GAMS/48 and solved using SCIP solver. To demonstrate engineering applicability, the thermodynamic modelling results are compared with Aspen EDR simulations, revealing relative errors ranging from 3.8 % to 15.6 %. The optimized design increases the overall heat transfer coefficient by 16.60 % and reduces the required area by 13.86 %. Furthermore, stability analysis under varying heat duties reveals clear operational boundaries, highlighting the importance of coordinating liquid level and skirt height to maintain feasible natural circulation.