Thermodynamic equilibria and concentrations in thermodynamic equilibria are of major importance in chemistry, chemical engineering, physical chemistry, medicine etc. due to a vast spectrum of applications. E.g., concentrations in thermodynamic equilibria play a central role for the estimation of drug delivery, the estimation of produced mass of products of chemical reactions, the estimation of deposited metal during electro plating and many more. Species concentrations in thermodynamic equilibrium are determined by the system of reactions and to the reactions’ associated stability constants. In many applications the stability constants and the system of reactions need to be determined. The usual way to determine the stability constants is to evaluate titration curves. In this context, many numerical methods exist. One major task in this context is that the corresponding inverse problems tend to be unstable, i.e., the output is strongly affected by measurement errors, and can output negative stability constants or negative species concentrations. In this work an alternative model for the species distributions in thermodynamic equilibrium, based on the models used for HySS or Hyperquad, and titration curves is presented, which includes the positivity of species concentrations and stability constants intrinsically. Additionally, in this paper a stabilized numerical methodology is presented to treat the corresponding model guaranteeing the convergence of the algorithm. The numerical scheme is validated with clinical numerical examples and the model is validated with a Citric acid–Nickel electrolyte. This paper finds a stable, convergent and efficient methodology to compute stability constants from potentiometric titration curves.