Bayesian Gaussian Process Optimization can be considered as a method for the determination of the model parameters, based on experimental data. In the range of soft QCD physics, the processes of hadron and nuclear interactions require using phenomenological models containing many parameters. In order to minimize the computation time, the model predictions can be parameterized using Gaussian Process regression, and then provide the input to the Bayesian Optimization. In this paper, the Bayesian Gaussian Process Optimization has been applied to the Monte Carlo model with string fusion. The parameters of the model are determined using experimental data on multiplicity and cross section of pp, pA and AA collisions in a wide energy range. The results provide important constraints on the transverse radius of the quark-gluon string ($r_{str}$) and the mean multiplicity per rapidity from one string ($\mu_0$).