A discrete adjoint-based method is employed to control multi-mode Rayleigh–Taylor (RT) instabilities via strategic manipulation of the initial interfacial perturbations. We seek to find to what extent mixing and growth can be enhanced at late stages of the instability and which modes are targeted to achieve this. Three objective functions are defined to quantify RT mixing and growth: (i) variance of mole fraction, (ii) a kinetic energy norm based on the vertical velocity and (iii) variations of mole fraction with respect to the unperturbed initial state. The sensitivity of these objective functions to individual amplitudes of the initial perturbations are studied at various stages of the RT instability. The most sensitive wavenumber during the early stages of the instability closely matches the most unstable wavenumber predicted by linear stability theory. It is also shown that randomly changing the initial perturbations has little effect at early stages, but results in large variations in both RT growth and its sensitivity at later times. The sensitivity obtained from the adjoint solution is employed in gradient-based optimization to both suppress and enhance the objective functions. The adjoint-based optimization framework was capable of completely suppressing mixing by shifting all of the interfacial perturbation energy to the highest modes such that diffusion dominates. The optimal initial perturbations for enhancing the objective functions were found to have a broadband spectrum resulting in non-trivial coupling between modes and depends on the particular objective function being optimized. The objective functions were increased by as much as a factor of nine in the self-similar late-stage growth regime compared to an interface with a uniform distribution of modes, corresponding to a 32% increase in the bubble growth parameter and 54% increase in the mixing width. It was also found that the interfacial perturbations optimized at early stages of the instability are unable to predict enhanced mixing at later times, and thus optimizing late-time multi-mode RT instabilities requires late-time sensitivity information. Finally, it was found that the optimized distribution of interfacial perturbations obtained from two-dimensional simulations was capable of enhancing the objective functions in three-dimensional flows. As much as 51% and 99% enhancement in the bubble growth parameter and mixing width, respectively, was achieved, even greater than what was reached in two dimensions.