Current multi-objective evolutionary polynomial regression (EPR) methodology has difficulties on decision-making of optimal EPR model. This paper proposes an intelligent multi-objective optimization-based EPR technique with multi-step automatic model selection procedure. A newly developed multi-objective differential evolution algorithm (MODE) is adopted to improve the optimization performance. The proposed EPR process is composed of two stages: (1) intelligent roughing model selection and (2) model delicacy identification. In the first stage, besides two objectives (model accuracy and model complexity), the model robustness measured by robustness ratio is considered as an additional objective in the multi-objective optimization. In the second stage, a new indicator named selection index is proposed and incorporated to find the optimal model. After intelligent roughing selection and delicacy identification, the optimal EPR model is obtained considering the combined effects of correlation coefficient, size of polynomial terms, number of involved variables, robustness ratio and monotonicity. To show the practicality of the proposed EPR technique, three illustrative cases helpful for geotechnical design are presented: (a) modelling of compressibility, (b) modelling of undrained shear strength and (c) modelling of hydraulic conductivity. For each case, a practical formula with better performance in comparison with various existing empirical equations is finally provided. All results demonstrate that the proposed intelligent MODE-based EPR technique is efficient and effective.