In water pollution control of urban river basins, one major issue is to reduce water pollutants at the least costs. In this research, an integrated hydrodynamic-water quality simulation-based multi-response surface method was proposed to help optimize mitigation measures on water quality rehabilitation in Shiwuli River basin of Hefei City, China. Initially, multi-response surface method (MRSM) with two-levels and two responses, i.e., cost-effectiveness ratio (CER) and water quality standard-reaching rate (WQSR) were designed where the factors of interests include the point source reducing rate (i.e., Factor A), the non-point source reducing rate (i.e., Factor B) and water transfer (i.e., Factor C). Then a hydrodynamic-water quality model was developed to predict water quality effects of the proposed combined factors. At the same time, the cost functions were established to estimate economic costs of the combined factors. Lastly, the analysis of variance (ANOVA) and multiple regression analysis (MRA) were employed to help gain some new cognition in the statistically significant way. The obtained model revealed that Factor A and B were significant for both CER and WQSR. This means that in the study area, pollution source reduction is more important than water transfer. Though the optimal minimal CER was obtainable with 0.1493 (M$ per percentage) when the pointing source reducing rate is 77.9%, the non-point source reducing rate is 40% and the water transfer was controlled at a low level, respectively, with water quality standard-reaching rate 70.7%. It fails to achieve water quality objectives (WQSR ≥ 90.0%). Moreover, In order to meet water quality target, simultaneous optimization of two objective functions: minimal CER and maximal WQSR was achievable with the optimum value of 90% for WQSR and 0.1810 (M$ per percentage) for CER ,respectively, with the pointing source reducing rate (90.0%), the non-point source reducing rate (57.2%) and the water transfer at a low level. It provided an alternative simulation–optimization method for decision of water pollution control which had the advantages in reflecting the real water environment system and simplifying the algorithm compared to the conventional mathematical programming at the watershed scale.