Abstract: N-sided hole filling plays an important role in surface modeling. To address twist compatibility issues, an algorithm based on constrained least square method to fill N-sided holes with NURBS patches is proposed. Firstly, initial cross-boundary derivatives of the surface are generated through linear interpolation of boundary derivatives. Then, a central point is calculated, and this central point along with the midpoints of the boundary lines are used as the first and last control points to construct internal lines, incorporating tangent vectors at the endpoints. Subsequently, normal constraints and twist compatibility conditions are introduced, and the linear interpolation energy of the initial surface boundary derivative are selected as the optimization objective function to generate the cross-boundary derivatives that satisfy boundary compatibility, twist compatibility and ε-G¹ continuity, where ε is the angular deviation given by the user. Finally, bicubic blended Coons surfaces are employed to achieve G¹ continuous filling for each rectangular region of the N-sided hole. The proposed method solves the numerical instability and surface smoothing problems caused by the small knot vector interval introduced by the Piegl’s algorithm to meet the twist compatibility. The angular deviation between normal vectors of adjacent surfaces is controlled within the user-specified ε angle tolerance. To validate the algorithm’s effectiveness, five representative experimental datasets, each emphasizing distinct aspects of size and curvature, were selected.