To accurately determine phase boundaries and phase transitions, thermodynamic models that describe phase free energies often have to be optimized based on experimentally observed phase equilibria. While different approaches exist for thermodynamic optimization, they are often implemented in ways that are not compatible with machine learning workflows requiring differentiable calculation of the loss function. In this work, we derive a phase-equilibrium loss function based on thermodynamic potentials that can be efficiently evaluated and enables gradient-based optimization by auto-differentiation in the PyTorch package. By minimizing this loss function, general thermodynamic model parameters can be optimized with respect to experimental phase-equilibria data. Using thermodynamic models in the CALculation of PHAse Diagram (CALPHAD) framework, we illustrate successful and efficient optimization in different systems, including ternary systems with more than 100 parameters. As the loss function is defined independently of the details of the thermodynamic models, it can be used to optimize machine learning thermodynamic models in general. In particular, we demonstrate top-down optimization of an atomistic potential from target phase equilibria.