In this paper we describe recent developments in the application of mathematical and computational techniques to the problem of designing binary gratings on top of a multilayer stack in such a way that the propagating modes have a specified intensity or phase pattern for a chosen range of wavelengths or incidence angles. The diffraction problems are transformed to strongly elliptic variational formulations of quasi periodic transmission problems for the Helmholtz equation in a bounded domain coupled with boundary integral representations in the exterior. We obtain analytic formulae for the gradients of cost functionals with respect to the parameters of the grating profile and the thickness of the layers, so that the optimal design problems can be solved by minimization algorithms based on gradient descent. For the computation of diffraction efficiencies and gradients the variational problems are solved by using a generalized finite element method with minimal pollution. We provide some numerical examples to demonstrate the convergence properties for evaluating diffraction efficiencies and gradients. The method is applied to optimal design problems for polarisation gratings and beam splitters.