We study the numerical approximation of convex optimal control problems governed by elliptic partial differential equations with oscillating coefficients. Since the objective functional contains flux, we approximate the problems using the mixed finite element methods. We first analyze the standard finite element approximation schemes. Then, motivated by the numerical simulation of the primal variable and the flux in highly heterogeneous porous media, we use a mixed multiscale finite element method for solving the state equations. The multiscale finite element bases are constructed by locally solving Dirichlet boundary value problems. The analysis of the approximate control problems is carried out under the assumption that the oscillating coefficients are locally periodic, which allows us to use homogenization theory to obtain the asymptotic structure of the solutions, although the numerical schemes are designed for general cases.