3D printing technology has been developing rapidly in recent years, but product quality control has become one of the main obstacles to its widespread use in manufacturing. A new stochastic computer model and robust optimization method are proposed for the highly fluctuating 3D printing process to improve the stability of the printed product quality. Firstly, the signal and noise are jointly modeled, and the idea of latent variables in machine learning is incorporated to overcome the limitation that the replication times of the stochastic Kriging model must be greater than one. Then, the chain rule and Woodbury identity are utilized to reduce the time required for hyperparameter estimation of the model. Finally, the optimization objective function is constructed based on the Taguchi quality loss function, and optimal process parameters are found using a genetic algorithm. The numerical simulation results demonstrate that the robust optimization method based on heteroskedasticity Gaussian process model proposed in this paper can estimate model hyperparameters faster and predict results more accurately. Furthermore, the prediction and validation results of 3D printing experiments verify the effectiveness of the proposed method.