This paper addresses the online modeling for time-varying manufacturing systems with random unknown model variations between production batches. By modeling the system as a Gaussian process, we first apply the standard Gaussian process regression (GPR) method for estimating the system model, which provides the optimal model estimate with the minimum mean square error (MSE). Then, an iterative form of the method is derived which is more computation efficient but maintains the estimation optimality. However, such optimality is obtained by continuously updating the covariances between the estimated model values and the measurements, which would make the storage and computation unaffordable when the control input can vary within an infinite control space. Due to such a limitation, a suboptimal interactive GPR method is further proposed by trading off the computation efficiency and the estimation accuracy, where the trade-off can be tuned by a designed parameter. Finally, effectiveness and performance of the proposed methods are demonstrated via both simulation and case study by comparing to the conventional nonlinear modeling methods. (C) 2016 Elsevier Ltd. All rights reserved.