Properly treating uncertainty is critical for robust system sizing of nearly/net zero energy buildings (ZEBs). To treat uncertainty, the conventional method conducts Monte Carlo simulations for thousands of possible design options, which inevitably leads to computation load that is heavy or even impossible to handle. In order to reduce the number of Monte Carlo simulations, this study proposes a response-surface-model-based system sizing method. The response surface models of design criteria (i.e., the annual energy match ratio, self-consumption ratio and initial investment) are established based on Monte Carlo simulations for 29 specific design points which are determined by Box-Behnken design. With the response surface models, the overall performances (i.e., the weighted performance of the design criteria) of all design options (i.e., sizing combinations of photovoltaic, wind turbine and electric storage) are evaluated, and the design option with the maximal overall performance is finally selected. Cases studies with 1331 design options have validated the proposed method for 10,000 randomly produced decision scenarios (i.e., users’ preferences to the design criteria). The results show that the established response surface models reasonably predict the design criteria with errors no greater than 3.5% at a cumulative probability of 95%. The proposed method reduces the number of Monte Carlos simulations by 97.8%, and robustly sorts out top 1.1% design options in expectation. With the largely reduced Monte Carlo simulations and high overall performance of the selected design option, the proposed method provides a practical and efficient means for system sizing of nearly/net ZEBs under uncertainty.