This article is devoted to developing an accurate operational matrix (OM) method for the solution of a new category of nonlinear optimal control problems (OCPs) explained by fractal‐fractional dynamical systems. The fractal‐fractional differentiation is considered in the sense of Atangana‐Riemann‐Liouville. The method is based on the shifted Legendre polynomials (LPs). Through the way, a new OM of fractal‐fractional differentiation is extracted for the shifted LPs. The optimality conditions are converted to an algebraic system of equations by using the shifted LP expansions of the state and control variables and applying the method of constrained extrema. More precisely, these variables are approximated by the shifted LPs with undetermined coefficients. Then, these expansions are replaced in the performance index, and the Gauss‐Legendre integration method is used to approximate the performance index for extracting a nonlinear algebraic equation. In the sequel, the mentioned OM and the collocation method are used to generate some algebraic equations from the dynamical system. Finally, the method of the constrained extrema is used by coupling the algebraic constraints generated from the dynamical system and the initial condition with the algebraic equation elicited from the performance index by a set of unknown Lagrange multipliers. The accuracy of the method is studied by solving some numerical examples.