<p>This paper addresses the numerical solution of delay fractional optimal control problems involving fractional derivatives in the Caputo sense. Fractional-order Chelyshkov functions (FCHFs) are employed to develop a spectral collocation method for solving these problems. The operational matrix of the Riemann–Liouville fractional integral for FCHFs is derived using the Laplace transform method and then used together with the Lagrange multiplier method to transform the fractional optimal control problem, which includes both state and control delays, into a system of algebraic equations in terms of the unknown FCHFs coefficients. An estimation of the error for the FCHFs approximation is provided and the convergence of the presented method is also discussed. The applicability and accuracy of the proposed method are demonstrated through various illustrative examples. Furthermore, some comparisons between our results and those published in the literature are presented. These numerical experiments show the superiority of the presented approach and that a small number of FCHFs is sufficient to achieve high-accuracy results.</p>