The authors consider an optimal control problem in the presence of uncertainty: the target function is the solution of an elliptic partial differential equation, steered by a control function and with a random field as input coefficient. They present a specially designed quasi-Monte Carlo method to approximate the expected values with respect to the uncertainty. Their method provides error bounds for the approximation of the stochastic integral, which do not depend on the number of uncertain variables. It results in faster convergence rates compared to Monte Carlo methods in the case of smooth integrands. The nonnegative quadrature weights preserve the convexity structure of the optimal control problem. The random field is, in principle, infinite-dimensional, and in practice, might need a large finite number of terms for accurate approximation. The novelty lies in using and analyzing a specially designed quasi-Monte Carlo method to approximate the possibly high-dimensional integrals with respect to the stochastic variables. The authors present a gradient descent algorithm to solve the optimal control problem for the case without control constraints and a projected variant of the algorithm for the problem with control constraints. Moreover, they present error estimates and convergence rates for the dimension truncation and the finite element discretization together with confirming numerical experiments.