[Objective] The multi-objective optimization of brushless direct current (BLDC) motors frequently relies on extensive finite element analysis, which is computationally expensive and time-consuming. While surrogate models accelerate the design process, traditional data-driven models often suffer from prediction distortion near physical boundaries, resulting in non-physical output values. To address this issue, this paper proposes a novel robust optimization methodology integrating a physically constrained gaussian process regression (GPR) model with a lower confidence bound (LCB) strategy. [Methods] Firstly, an improved Morris trajectory sampling method based on the Campolongo strategy was adopted to perform global sensitivity analysis on the key structural parameters of the motor, and achieved dimensionality reduction and screening of design variables. Secondly, considering the physical characteristics of motor efficiency and torque ripple, logarithmic (Log) and logistic (Logit) transformations were introduced into the GPR modeling to construct a high-accuracy surrogate model that satisfied physical boundary constraints. Finally, uncertainty quantification was performed using the predictive variance information provided by GPR, and a robust optimization strategy based on LCB was proposed. Through a two-stage search, the comprehensive performance of the motor was improved, while the engineering robustness of the design scheme was effectively enhanced. [Results] The proposed framework demonstrated excellent performance during validation. The improved GPR model mathematically eliminated non-physical predictions at the boundaries of motor efficiency and torque ripple, and it achieved significantly higher generalization accuracy on the test set compared with conventional surrogate models. When applied to the design of the BLDC motor, the final optimized design exhibited substantial improvements. Specifically, the output torque was increased by 9.56%, and the torque ripple was dramatically reduced by 47.18%. Meanwhile, the overall motor efficiency was reliably maintained above 90%. [Conclusion] This method effectively addresses the computational bottlenecks and model prediction distortions present in BLDC motor optimization. By structurally integrating physical constraints within the modeling architecture and utilizing the uncertainty-aware LCB strategy, the proposed method comprehensively improves electromagnetic performance while effectively enhancing the engineering robustness and practical viability of the design scheme.