The cooling load is the basis of cooling source design, while the inherent uncertainties of building information, weather condition, and internal load make it impossible to obtain a determinate load. Probability method can represent characteristics of uncertain cooling loads well, but a large number of Monte Carlo (MC) simulations should be carried out. The optimal cooling source design can be formulated as a mixed-integer-linear-programming model (MILP), which can be solved efficiently to obtain the global optimality using a state-of-the-art MILP solver. However, if all the MC simulations are used in the optimization problem, the size of MILP model would lead to computational issues or even be unsolvable. This study, therefore, proposes a robust optimization design model for sizing the cooling source when there is cooling load uncertainty. A method named in this paper cooling load bin is proposed and implemented by converting time series cooling loads obtained by MC simulations to those in the frequency domain. Therefore, the cooling load frequency and mean value in each cooling interval are obtained and used in the optimal design model, which can be solved efficiently by the General Algebraic Modeling System (GAMS). The robust model is applied to the optimal design of a cooling source to minimize the cost. A case study on the cooling source of a hospital in Tianjin is conducted to demonstrate the effectiveness of the proposed robust model. Furthermore, the accuracy of the solution and the computational efficiency used to evaluate the proposed robust model are systematically investigated and compared with the deterministic model.