This thesis is devoted to developing a robust Model Predictive Control (MPC) strategy based on Gaussian Processes (GP), especially for Drinking Water Networks (DWN). Nowadays there are many different MPC strategies developed for DWN, such as certain-equivalent MPC (CEMPC) and chance-constrained MPC (CC-MPC). The general control objectives for DWN are economic by managing the water supply to minimise water production and transport costs, all the tanks running in safe ways with their limitations and reducing the undesired abrupt changes by minimising their slew rate and obtaining smooth signals. For the deterministic system model, the control objectives are elementary fulfilled. But the main challenge for DWN is to propagate and incorporate exogenous and endogenous uncertainties to MPC closed loop over the prediction horizon. Considering the control-oriented model of the DWN, the MPC controller design is hereby divided into two parts: system disturbances forecasting and the robust MPC controller design. Case studies based on Barcelona DWN have been executed to verify the performance of proposed methodologies. The first part of this thesis leads to forecast system disturbances by using GP. In a DWN system, system disturbances come mainly water demands associated to consumer sectors. Hence, it is necessary to model each water demand and forecast the water demand in a short term that covers the MPC prediction horizon. GP regression is regarded as one of state-of-the-art regression methods able to select model parameters by using Bayesian Inference theory with a collection of past data. Besides, it is believed that the GP regression method has a difficult for the multiple-step ahead forecasting. Hence, the Double-seasonal Holt-winters method is used for forecasting the expected disturbances while the stochastic disturbances are forecasted by using GP. Finally, the desired forecasting results are a set of Gaussian distributions over the MPC prediction horizon. The second part of this thesis is to incorporate the forecasting results from GP within MPC closed loop. This MPC strategy based on GP is named GP-MPC. Using the given system model, the deterministic state evolutions can be obtained while the uncertainty of state propagation over a given prediction horizon can be also achieved though the linear approximation of GP. Therefore, the worst-case state evolutions over the MPC prediction horizon can also be determined in the MPC cost function and constraints. The desired performance of applying GP-MPC in the closed-loop system is that the system has more safety than the CE-MPC and meanwhile it probably brings more expenses. Comparisons of GP-MPC and previous developed approaches are carried out by a case study of the three-tank system inspired in the Barcelona DWN. A set of key performance indicators are defined to compare performances of different MPC strategies. Finally, through the simulation results, the GP-MPC has the similar performance as the CC-MPC, both of which have much more expenses than the CE-MPC. As a result of considering the uncertainties inside the system, more expenses is necessary to maintain the safety of the whole system. Hence, the GP-MPC is more advanced. Moreover, the proposed GP-MPC is required to be tested with the whole DWN and using the real data from a DWN system. So the future works of this thesis have been outlined.