A path to enhance the coefficient of performance (COP) of domestic heat pumps is to replace the throttling valve with a turboexpander. Designing this component requires addressing supersonic flashing flows localized within the nozzle. In particular, two-phase flows prevent the application of the method of characteristics to design efficient converging–diverging shapes. To overcome this issue, we propose a shape optimization combined with a homogeneous equilibrium flow model as a design tool. Kriging models are used as surrogates of the objective function and constraint to reduce the overall computational burden associated with shape optimization. The infilling criterion to improve surrogate accuracy involves the maximization of the constrained expected improvement. A multi-point optimization is formulated, including design and off-design conditions in both objective and constraint, for a total of four relevant operating conditions. The off-design conditions are derived from system analysis, accounting for variations in the evaporator temperature and outlet gas-cooler temperature. The optimized nozzle results in a COP improvement of 7.6%, compared to a 6.6% increase when replacing the throttling valve with a baseline geometry designed as recommended by the literature. The adoption of a multi-point optimization strategy ensures that this enhancement in design conditions is not offset by off-design performance decay, which is 0.3 percentage points over the three off-design conditions. Overall, replacing the throttling valve with the optimized expander leads to an average COP improvement of 8.6%.