  This manuscript presents a framework for using multilevel quadrature formulae
to compute the solution of optimal control problems constrained by random
partial differential equations. Our approach consists in solving a sequence of
optimal control problems discretized with different levels of accuracy of the
physical and probability discretizations. The final approximation of the
control is then obtained in a postprocessing step, by suitably combining the
adjoint variables computed on the different levels. We present a convergence
analysis for an unconstrained linear quadratic problem, and detail our
framework for the specific case of a Multilevel Monte Carlo quadrature formula.
Numerical experiments confirm the better computational complexity of our MLMC
approach compared to a standard Monte Carlo sample average approximation, even
beyond the theoretical assumptions.
