In this paper, a constrained distributed optimal control problem governed by a first-order elliptic system is considered. Least-squares mixed finite element methods, which are not subject to the Ladyzhenkaya-Babuska-Brezzi consistency condition, are used for solving the elliptic system with two unknown state variables. By adopting the Lagrange multiplier approach, continuous and discrete optimality systems including a primal state equation, an adjoint state equation, and a variational inequality for the optimal control are derived, respectively. Both the discrete state equation and discrete adjoint state equation yield a symmetric and positive definite linear algebraic system. Thus, the popular solvers such as preconditioned conjugate gradient (PCG) and algebraic multi-grid (AMG) can be used for rapid solution. Optimal a priori error estimates are obtained, respectively, for the control function in $L^2(Ω)$-norm, for the original state and adjoint state in $H^1(Ω)$-norm, and for the flux state and adjoint flux state in $H$(div; $Ω$)-norm. Finally, we use one numerical example to validate the theoretical findings.