Los equilibrios termodinámicos y las concentraciones en equilibrios termodinámicos son de gran importancia en química, ingeniería química, química física, medicina, etc. debido a un vasto espectro de aplicaciones. Por ejemplo, las concentraciones en equilibrios termodinámicos desempeñan un papel central para la estimación de la entrega de medicamentos, la estimación de la masa producida de productos de reacciones químicas, la estimación del metal depositado durante el electroplateado y muchos más. Las concentraciones de especies en equilibrio termodinámico son determinadas por el sistema de reacciones y por las constantes de estabilidad asociadas a las reacciones. En muchas aplicaciones, las constantes de estabilidad y el sistema de reacciones necesitan ser determinados. La forma habitual de determinar las constantes de estabilidad es evaluar las curvas de titulación. En este contexto, existen muchos métodos numéricos. Una tarea importante en este contexto es que los problemas inversos correspondientes tienden a ser inestables, es decir, la salida se ve fuertemente afectada por errores de medición, y puede producir constantes de estabilidad negativas o concentraciones de especies negativas. En este trabajo se presenta un modelo alternativo para las distribuciones de especies en equilibrio termodinámico, basado en los modelos utilizados para HySS o Hyperquad, y en curvas de titulación, que incluye la positividad de las concentraciones de especies y constantes de estabilidad de forma intrínseca. Además, en este artículo se presenta una metodología numérica estabilizada para tratar el modelo correspondiente garantizando la convergencia del algoritmo. El esquema numérico se valida con ejemplos numéricos clínicos y el modelo se valida con un electrolito de ácido cítrico-níquel. Este artículo encuentra una metodología estable, convergente y eficiente para calcular constantes de estabilidad a partir de curvas de titulación potenciométrica.