Abstract In this study, a copula-based flexible-stochastic programming (CFSP) method is developed for planning regional energy system (RES). CFSP can deal with multiple uncertainties expressed as interval values, random variables and fuzzy sets as well as their combinations employed to objective function and soft constraints. It can also reflect uncertain interactions among random variables through using copula functions even having different probability distributions and previously unknown correlations. Then, based on the developed CFSP approach, a CFSP-RES model is formulated for planning RES of the urban agglomeration of Beijing and Tianjin (China). Results disclose that uncertainties existed in the system components have significant effects on the outputs of decision variables and system cost, and the variation of system cost is reached 16.3%. Results also reveal that air pollutant emissions can be mitigated if the urban agglomeration can co-implement renewable energy development plans (REDP) over the planning horizon, with the reductive rates of [3.3, 7.6] % of sulfur dioxide (SO2), [2.7, 4.1] % of nitrogen oxides (NOx) and [7.0, 11.5] % of particulate matter (PM10). Compared to joint-probabilistic chance-constrained programming (JCP), the CFSP method is more effective for handling multiple random parameters associated with different probability distributions in which their correlations are unknown. Thus, it is not limited to some unjustified assumptions and can be applied to a wider range of problems than previous studies. The findings are helpful to explore the influence of interaction among random variables on modeling outputs and provide in-depth analysis for identifying desired decision schemes for planning RES.