enter image description hereI’m having an issue using the least_squares function to solve a system of non-linear equation to obtain a chemical speciation (here of the system Nickel-ammonia). The system of equations is given here :

def equations(x):
    Ni0, Ni1, Ni2, Ni3, Ni4, Ni5, Ni6, NH3 = x
    eq1 = np.log10(Ni1)-np.log10(Ni0)-np.log10(NH3)-K1
    eq2 = np.log10(Ni2)-np.log10(Ni1)-np.log10(NH3)-K2
    eq3 = np.log10(Ni3)-np.log10(Ni2)-np.log10(NH3)-K3
    eq4 = np.log10(Ni4)-np.log10(Ni3)-np.log10(NH3)-K4
    eq5 = np.log10(Ni5)-np.log10(Ni4)-np.log10(NH3)-K5
    eq6 = np.log10(Ni6)-np.log10(Ni5)-np.log10(NH3)-K6
    eq7 = Ni0 + Ni1 + Ni2 + Ni3 + Ni4 + Ni5 + Ni6 - Ni_tot
    eq8 = Ni1 + 2*Ni2 + 3*Ni3 + 4*Ni4 + 5*Ni5 + 6*Ni6 - NH3_tot[i]
    
    return np.array([eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8])

I use a loop to calculate the values of the parameters at various NH3_tot, allowing to observe the variation of concentrations of the various species. The solving is done using the following :

# Arrays for data storage
NH3_tot = np.arange(0, 6, 0.1)
Ni0_values = []
Ni1_values = []
Ni2_values = []
Ni3_values = []
Ni4_values = []
Ni5_values = []
Ni6_values = []
Solutions = []

# Speciation determination with least_square at various NH3_tot
for i in range(len(NH3_tot)):
    x0 = [1,0.1,0.1,0.1,0.1,0.1,0.1,0.1]
    result = least_squares(equations, x0, bounds=([0,0,0,0,0,0,0,0],[Ni_tot,Ni_tot,Ni_tot,Ni_tot,Ni_tot,Ni_tot,Ni_tot,Ni_tot]))
    
    
    Solution = result.x
    Solutions.append(Solution)
    Ni0_value,Ni1_value,Ni2_value,Ni3_value,Ni4_value,Ni5_value,Ni6_value, NH3 = Solution
    
    # Extract concentrations of species
    Ni0_values.append(Ni0_value)[enter image description here](https://i.sstatic.net/pQaZ1efg.png)
    Ni1_values.append(Ni1_value)
    Ni2_values.append(Ni2_value)
    Ni3_values.append(Ni3_value)
    Ni4_values.append(Ni4_value)
    Ni5_values.append(Ni5_value)
    Ni6_values.append(Ni6_value)

The issue I get is that for some values (generally the ones for low and high NH3_tot values), are completelly wrong, especially not respecting eq. 7, as the sum of the species exceeds Ni_tot.

I tried moving the bounds, but it does not seem to help, and the extrem values are generally wrong.

Printing the result always shows :
message: xtol termination condition is satisfied.
success: True

Indicating that the function effectively operates correctly.

Changing the chemical system results in similar output.

Changing the writting of the equations impacts greatly the result (If know why please help !!), and the use of np.log10 gives the most accurate, but still unsatisfactory.

Changing the method of resolution of the least_squares function does not improve the result.