The performance of supersonic ejectors for heat recovery refrigeration cycles is optimised by means of fluid dynamic shape optimisation with a discrete adjoint method to efficiently evaluate gradients of an objective function with respect to an arbitrary number of design variables. The Reynolds-averaged Navier-Stokes equations and the discrete adjoint equations are solved using the finite volume solver SU2. The model is validated on cases from the literature. It is conclusively demonstrated that the method employed can generate a well performing ejector shape from a failed design despite concerns about the well-posedness of the problem given the apparent discontinuity of the entrainment ratio objective function at the critical point in the parameter space. This is achieved in only 10-20 iterations with a cost of one or two direct solutions each. The entrainment ratios predicted with the optimal shape exceed those announced by a leading manufacturer for identical cycle conditions by around 15 %. More importantly, discrete adjoint shape optimisation is here applied successfully for the first time to transonic ejector flows and is shown to be reliable and robust, paving the way for the application of topology optimisation methods to ejectors.