We employ the Galerkin Mittag-Leffler method to address multi-dimensional fractional optimal control problems (MFOCPs) in the sense of the Caputo fractional derivative, subject to both equality and inequality constraints. The Riemann–Liouville operational matrix for Mittag-Leffler functions is derived and utilized alongside the Galerkin method to express the MFOCP as a system of algebraic equations that facilitates efficient computation. The properties of convergence and the error estimation for the Mittag-Leffler polynomials are thoroughly examined. Additionally, the proposed method is tested on four examples, and the results are compared to those reported by other researchers. The comparisons confirmed the accuracy and effectiveness of our approach. In certain cases, the solutions obtained were exact, demonstrating the robustness and precision of the method.