In this paper the method of penalized high order Legendre polynomials and specific designs of experiments for model building on the basis of numerical or physical experiments is proposed. The generalized thin-plate energy functional is used for penalization of the least square functional. The optimal choice of smoothing parameter (penalization coefficient) is implemented using the cross-validation method. Special axial-symmetrical and 90 degree rotational symmetrical D-optimal experimental designs are obtained using the direct constrained optimization method. The method was tested for known optimization test problems with 2–5 variables and showed prediction accuracy comparable with kriging.