Numerical solver finds different answer depending on guess

Hi,

I made this brief code for determining the distribution of different molecules at a given pH. The problem is, that the solver is extremely dependent on the initial guess I provide it. I have tried with a variety of solver, but could not figure out how to optimally determine the solution to the four equations. Previously I would solve these equations symbolic in e.g. Maple, but I do not have that tool available at my current work. I know I could just isolate the expressions, but this is also relevant for more complicated expressions. How can I find a solution to this function? The pH value can vary between 0 and 14, where it is much more difficult to solve between e.g. 6 and 10.

I have access to matlab, control system, optimization and system identification

%% Calculate the number of protons a single soda molecule binds with. The solver may not find a solution if guess is way off

p.k_H = 29.76;                                                              % Henry constant
p.K_a1 = 2.5*10^-4;                                                         % Equilibrium constant 1 (H2CO3 <-> HCO3- + H+)
p.K_a2 = 4.69*10^-11;                                                       % Equilibrium constant 2 (HCO3- <-> CO3-- + H+)
p.K_h = 10^-2.8;                                                            % Hydration constant (H2CO3/CO2)

pH = 4.92;                                                                  % pH-value

r0 = [10^-5 10^-6 10^-7 10^-8]/10;                                         % Initial guess

options = optimset('TolFun',10^-20,'MaxFunEvals',1000,'MaxIter',1000);    % Solver options
val = lsqnonlin(@(r) myfun(r,pH,p),r0,zeros(4,1),[],options)                % Molar concentrations

eq = sum(val.*[2 2 1 0])/(sum(val));                                        % How many protons a Na2CO3 molecule absorbs at this pH

fprintf('Na2CO3 is equivalent to %.3f molar units of base\n',eq)

function F = myfun(r,pH,p)

c_CO2 = r(1);
c_H2CO3 = r(2);
c_HCO3 = r(3);
c_CO3 = r(4);

c_H = 10^-pH;

F(1) = 10^-14/c_H + c_HCO3 + 2*c_CO3 - c_H;                                 % Charge neutrality
F(2) = c_H2CO3/c_CO2 - p.K_h;                                               % Hydration
F(3) = c_H*c_CO3/c_HCO3 - p.K_a2;                                           % Eq2
F(4) = c_H*c_HCO3/c_H2CO3 - p.K_a1;                                         % Eq1

end

Answer 1

[John D'Errico]

This question suffers from several misunderstandings about optimization tools, several common problems.

First, think of an optimizer as a blindman, placed somewhere on the face of the earth. Give him only his cane, then task him with finding the lowest spot on earth. Yes, I hope you will give him scuba gear too, as that spot will be deep in a Pacific ocean trench. But the point is, if you set him down to start with in the Himalayas, do you exeect him to find the Marianas trench in the Pacific? As bad, supose you set him down near the Dead Sea? He will find a low spot, but then decide that no direction leads him down farther, so he will stop searching.

An optimization tool is the same. If you start it in the wrong place, then what do you expect? It looks aroun, then may get lost if you give it different starting values. All is well and good when a problem has only ONE local minimizer, so that solution is also the globally best solution. But many problems do not have this property. So change the starting values, and what do you expect? Good starting values often matter.

Next, problems that vary over many powers of 10 in the paameters are nasty things. Worse, your problem goes to hell if the optimizer looks even slightly below zero, yet your parameters are very near zero in absolute terms.

The good thing is this is easily repaired. Work in terms of log to the base 10. So define the unknowns as log10 of your unknowns.

Your initial guess now becomes

r0 = [-5 -6 -7 -9] - 1;
% Initial guess

Inside your objective, the very first thing you will do is raise 10 to the power of r.

function F = myfun(r,pH,p)
R = 10.^r;
c_CO2 = R(1);
c_H2CO3 = R(2);
c_HCO3 = R(3);
c_CO3 = R(4);
...

But as far as the optimizer is concerned, now all the unknowns have the same order of magnitude, and things are now well scaled. As well, it never has to worry about numbers going less than zero, since the transformation handles everything.

Oh, and trying to set TolFun to 1e-20? Silly. You can't achieve that much accuracy in a problem like this.

So, first, I'll hack your objective function a bit, just to avoid needing to pass in those parameters. I am feeling lazy today.

function F = myfun(r)
  p.k_H = 29.76;         % Henry constant
  p.K_a1 = 2.5*10^-4;      % Equilibrium constant 1 (H2CO3 <-> HCO3- + H+)
  p.K_a2 = 4.69*10^-11;    % Equilibrium constant 2 (HCO3- <-> CO3-- + H+)
  p.K_h = 10^-2.8;         % Hydration constant (H2CO3/CO2)
  pH = 4.92;               % pH-value
  
  R = 10.^r;
  
  c_CO2 = R(1);
  c_H2CO3 = R(2);
  c_HCO3 = R(3);
  c_CO3 = R(4);
  c_H = 10^-pH;
  
  F(1) = 10^-14/c_H + c_HCO3 + 2*c_CO3 - c_H;  % Charge neutrality
  F(2) = c_H2CO3/c_CO2 - p.K_h;                % Hydration
  F(3) = c_H*c_CO3/c_HCO3 - p.K_a2;            % Eq2
  F(4) = c_H*c_HCO3/c_H2CO3 - p.K_a1;          % Eq1

Next, test it!!!!!!!! Verify that it does something reasonable.

myfun([-6 -7 -8 -9])
ans =
   -1.201e-05     0.098415   1.2022e-06   -0.0002488

So it produces numbers, 4 of them as expected, varying over 4 orders of magnitude. That in itself will be highly problematic. Why? A tool like lsqnonlin tries to minimize the sum of squares of the numbers it sees. So if one of those values is 4 powers of 10 larger than the others, it can essentially ignore the other three. All that matters is the biggest one.

So merely using lsqnonlin here MAY be a highly suspect thing to do. And settign TolFun at 1e-20? I'd bet it is way too optimistic. But for that tight a tolerance, you will need to increase the maximum function evals allowed. (And learn to use scientific notation with your numbers.)

Next, set the display parameter to iter, so we can see if it is converging. ALWAYS do this the first time you run a specific optimization.

options = optimset('TolFun',1e-15,'MaxFunEvals',5000,'MaxIter',1000,'Display','iter');

We can now try lsqnonlin, but there is no longer a need for a lower bound, because we are working in a log space instead.

r0 = [-6 -7 -8 -9];
val = lsqnonlin(@myfun,r0,[],[],options)

                                         Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality
     0          5       0.0096856                        0.0227
     1         10      0.00129798             10        0.00308      
     2         15      0.00129798        5.31272        0.00308      
     3         20     0.000173209        1.32818       0.000421      
     4         25     2.00996e-05       0.948856       5.85e-05      
     5         30     1.78119e-06        1.02082       8.18e-06      
     6         35     8.57658e-08        1.05491       1.01e-06      
     7         40     3.71167e-09        1.06807       5.56e-08      
     8         45     3.59717e-10        1.00663       1.35e-09      
     9         50     3.54168e-10     0.00409291       6.69e-10      
    10         55     3.48651e-10     0.00591474       1.33e-09      
    11         60     3.43293e-10     0.00399236        6.4e-10      
    12         65     3.37891e-10     0.00597722       1.31e-09      
    13         70     3.32645e-10     0.00389634       6.26e-10      
    14         75      2.0284e-11       0.828326       1.77e-10      
    15         80     2.01912e-11    0.000524991       1.06e-10      
    16         85     2.00988e-11    0.000775195       1.76e-10      
    17         90     2.00068e-11     0.00052201       1.05e-10      
    18         95     1.99152e-11    0.000772959       1.75e-10      
    19        100     1.98242e-11    0.000519053       1.04e-10      
    20        105     1.97335e-11    0.000770741       1.75e-10      
    21        110     1.96433e-11    0.000516122       1.04e-10      
    22        115     1.95534e-11    0.000768542       1.74e-10    
    
    ...
        
   980       4905     1.12224e-14    8.67351e-05       7.73e-12      
   981       4910     1.10723e-14    1.77985e-05       1.27e-12      
   982       4915     1.09241e-14    8.56892e-05       7.63e-12      
   983       4920     1.07778e-14    1.75701e-05       1.25e-12      
   984       4925     1.06333e-14    8.46533e-05       7.53e-12      
   985       4930     1.04908e-14    1.73441e-05       1.24e-12      
   986       4935       1.035e-14    8.36273e-05       7.43e-12      
   987       4940     1.02111e-14    1.71207e-05       1.22e-12      
   988       4945     1.00739e-14    8.26112e-05       7.34e-12      
   989       4950     9.93853e-15    1.68998e-05        1.2e-12      
   990       4955     9.80485e-15     8.1605e-05       7.24e-12      
   991       4960     9.67296e-15    1.66814e-05       1.19e-12      
   992       4965     9.54272e-15    8.06085e-05       7.15e-12      
   993       4970     9.41422e-15    1.64655e-05       1.17e-12      
   994       4975     9.28731e-15    7.96218e-05       7.06e-12      
   995       4980     9.16212e-15     1.6252e-05       1.15e-12      
   996       4985     9.03848e-15    7.86448e-05       6.96e-12      
   997       4990     8.91651e-15    1.60409e-05       1.14e-12      
   998       4995     8.79606e-15    7.76776e-05       6.87e-12      
   999       5000     8.67723e-15    1.58322e-05       1.12e-12      
  1000       5005     8.55989e-15      7.672e-05       6.78e-12  

  val =
      -3.4349      -6.2349      -4.9169       -8.481

I skipped the middle 900 or so iterations in that list.

What did that get us to?

myfun(val)
ans =
   2.9286e-06    2.623e-09   2.6118e-09   4.1699e-08

Trying to set the tolerance even that low is surely twiddling with how many angels can fit on the head of a pin. But lsqnonlin seems to work reasonably here. The biggest improvement was surely moving to the logs of your parameters, as that made the problem well scaled in the eyes of lsqnonlin.