Curve and surface fitting are important and attractive problems in many applied domains, from CAD techniques to geological prospections. Different methodologies have been developed to find a curve or a surface that best describes some 2D or 3D data, or just to approximate some function of one or several variables. In this paper, a new methodology is presented for optimal knots’ placement when approximating functions of one or two variables. When approximating, or fitting, a surface to a given data set inside a rectangle using B-splines, the main idea is to use an appropriate multi-objective genetic algorithm to optimize both the number of random knots and their optimal placement both in the x and y intervals, defining the corresponding rectangle. In any case, we will use cubic B-splines in one variable and a tensor product procedure to construct the corresponding bicubic B-spline basis functions in two variables. The proposed methodology has been tested both for functions of one or two independent variables, in order to evaluate the performance and possible issues of the procedure.