Abstract: N-sided hole filling plays a critical role in surface modeling and corner transition. To address twist compatibility issues, an algorithm based on constrained least squares energy optimization to fill N-sided holes with NURBS patches is proposed, based on Piegl’s algorithm. Initial cross-boundary derivatives of the surface are generated through linear interpolation of boundary derivatives. Normal constraints and twist compatibility conditions are introduced, and the smoothing energy of the cross-boundary derivative curve and the linear interpolation energy of the initial surface boundary derivative are selected as the optimization objective function to generate surface cross-boundary derivatives that satisfy boundary compatibility, twist compatibility and ε-G¹ continuity, where ε is the angular deviation given by the user. The proposed method solves the numerical instability and surface smoothing problems caused by the small node vector interval introduced by the Piegl’s algorithm to meet the twist compatibility. The algorithm can generate high-precision cross-boundary derivative curves that meet the requirements, and the angular error between the cross-boundary derivatives and the normal is within the ε angle tolerance range defined by the user.