Python - Nelder Mead vs fsolve: pH equilibria optimization problem leads to numerical instability in Nelder Mead?

I'm creating a pH equilibrium calculator - starting with initial compositions of certain acids / bases the calculator should come up with the equilibrated concentrations of each component based on physical constants.

I've formulated this as a set of equations which need to be solved simultaneously.

One key constraint is that H_plus*OH_minus must be 10^-14, and I've included this in the equations.

Code #1: Fsolve. This works, but is sensitive to initial conditions, and I can't put constraints (I need concentrations to be >=0).

import numpy as np
from scipy.optimize import fsolve

def equations(variables, Ka, Kw, initial_concentrations):
    HA, A_minus, H_plus, OH_minus, ion_Na = variables
    HA_init, A_minus_init, H_plus_init, OH_minus_init, ion_Na_init = initial_concentrations
    # Define the equations
    eq1 = HA + A_minus - HA_init - A_minus_init # Mass balance on acetate
    eq2 = H_plus + ion_Na - OH_minus - A_minus  # charge balance
    eq3 = H_plus * OH_minus - Kw                # Kw for water
    eq4 = Ka * HA - A_minus * H_plus            # HA acid dissociation
    eq5 = ion_Na - ion_Na_init                  # inert ion
    return [eq1, eq2, eq3, eq4, eq5]

def equilibrate(initial_concentrations, initial_guess):
    # User-supplied equilibrium constant (Ka) and initial concentration of acetic acid (HA_initial)
    Ka = 1.8e-5             # mol/L, acetic acid dissociation constant
    Kw = 1e-14              # mol/L, ionic product of water
    
    # Use fsolve from SciPy to find the equilibrium concentrations
    solution = fsolve(equations, initial_guess, args=(Ka, Kw, initial_concentrations), xtol=1e-18)
    
    return solution


initial_concentrations = [0.0, 0.1, 1e-7, 1e-7, 0.1]

initial_guess = initial_concentrations

# Extract the concentrations
HA, A_minus, H_plus, OH_minus, ion_Na = equilibrate(initial_concentrations, initial_guess)

# Calculate pH
pH = -np.log10(H_plus)
pKw = -np.log10(H_plus) -np.log10(OH_minus)

# Display the results
print(f"HA: {HA:.8f}")
print(f"A-: {A_minus:.8f}")
print(f"H+: {H_plus:.8f}")
print(f"OH-: {OH_minus:.8f}")
print(f"pH: {pH:.2f}")
print(f"pKw: {pKw:.2f}") # THIS IS THE CONDIRMATORY CHECK - Should be 14.

Code #2: Nelder Mead. This converges, despite the tight tolerances I set, but gives meaningless values - pKw

import numpy as np
from scipy.optimize import minimize

def equations(variables, Ka, Kw, initial_concentrations):
    HA, A_minus, H_plus, OH_minus, ion_Na = variables
    HA_init, A_minus_init, H_plus_init, OH_minus_init, ion_Na_init = initial_concentrations
    # Define the equations
    eq1 = HA + A_minus - HA_init - A_minus_init # Mass balance on acetate
    eq2 = H_plus + ion_Na - OH_minus - A_minus  # Charge balance
    eq3 = H_plus * OH_minus - Kw                # Kw for water
    eq4 = Ka * HA - A_minus * H_plus            # HA acid dissociation
    eq5 = ion_Na - ion_Na_init                  # Inert ion
    return eq1**2 + eq2**2 + eq3**2 + eq4**2 + eq5**2  # Use squared error for the minimization

def equilibrate(initial_concentrations, initial_guess, Ka, Kw):
    # Define bounds for variables (concentrations should be >= 0)
    bounds = [[0, 1], [0, 1], [0, 1], [0, 1], [0, 1]]
    
    # Use minimize from SciPy with BFGS method to find the equilibrium concentrations
    result = minimize(equations, initial_guess, args=(Ka, Kw, initial_concentrations),
                      bounds=bounds,method='Nelder-Mead', options={'disp': True, 'xatol':1E-40, 'fatol':1E-40, 'maxfev':1000000})
    
    return result.x

# User-supplied equilibrium constant (Ka) and initial concentration of acetic acid (HA_initial)
Ka = 1.8e-5             # mol/L, acetic acid dissociation constant
Kw = 1e-14              # mol/L, ionic product of water

initial_concentrations = [0.00, 0.1, 1e-7, 1e-7, 0.1]

initial_guess = initial_concentrations

# Extract the concentrations
HA, A_minus, H_plus, OH_minus, ion_Na = equilibrate(initial_concentrations, initial_guess, Ka, Kw)

# Calculate pH
pH = -np.log10(H_plus)

pKw = -np.log10(H_plus) -np.log10(OH_minus)

# Display the results
print(f"HA: {HA:.8f}")
print(f"A-: {A_minus:.8f}")
print(f"H+: {H_plus:.8f}")
print(f"OH-: {OH_minus:.8f}")
print(f"pH: {pH:.2f}")
print(f"pKw: {pKw:.2f}") # THIS IS THE CONDIRMATORY CHECK - Should be 14.

The confirmatory check value is all over the place, instead of being 14.

So my questions are:

1. Why does the Nelder Mead algorithm do so poorly? Is it because the absolute values involved are close to the error? I thought I could resolve this by tightening tolerances

2. Is there a way to impose positive constraints in fsolve?

3. Is there generally a better way to do this?

source https://stackoverflow.com/questions/77205188/python-nelder-mead-vs-fsolve-ph-equilibria-optimization-problem-leads-to-nume