Solving optimization problems with partial differential equations constraints is one of the most challenging problems in the context of industrial applications. In this paper, we study the finite volume element method for solving the elliptic Neumann boundary control problems. The variational discretization approach is used to deal with the control. Numerical results demonstrate that the proposed method for control is second-order accuracy in the <em>L</em><sup>2</sup> (Γ) and <em>L</em><sup>∞</sup> (Γ) norm. For state and adjoint state, optimal convergence order in the <em>L</em><sup>2</sup> (Ω) and <em>H</em><sup>1</sup> (Ω) can also be obtained.