This paper presents a novel tool for the shape optimization of turbomachinery blade profiles operating with fluids in non-ideal thermodynamic conditions and in complex flow configurations. In novel energy conversion systems, such as organic Rankine cycles or supercritical CO2 cycles, the non-conventional turbomachinery layout as well as the complex thermodynamics of the working fluid complicate significantly the blade aerodynamic design. For such applications, the design of turbomachinery may considerably benefit from the use of systematic optimization methods, especially in combination with high-fidelity computational fluid dynamics (CFD), as it is shown in this paper. The proposed technique is implemented in the shape-optimization package FORMA (Fluid-dynamic OptimizeR for turbo-Machinery Aerofoils) developed in-house at the Politecnico di Milano. FORMA is constructed as a combination of a generalized geometrical parametrization technique based on B-splines, a CFD solver featuring turbulence models and arbitrary equations of state, and multiple surrogate-based evolutionary strategies based on either trust-region or training methods. The application to the re-design of a supersonic turbine nozzle shows the capabilities of applying a high-fidelity optimization, consisting of a 50% reduction in the cascade loss coefficient and in an increased flow uniformity at the inlet of the subsequent rotor. Two alternative surrogate-based evolutionary strategies and different fitness functions are tested and discussed, including nonlinear constraints within the design process. The optimization study reveals relevant insights into the design of supersonic turbine nozzles as well on the performance, reliability, and potential of the proposed design technique.