Numerical solver finds different answer depending on guess

4 ビュー (過去 30 日間)
Bastian Andersen
Bastian Andersen 2018 年 12 月 4 日
編集済み: John D'Errico 2018 年 12 月 4 日
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

採用された回答

John D'Errico
John D'Errico 2018 年 12 月 4 日
編集済み: John D'Errico 2018 年 12 月 4 日
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.

その他の回答 (0 件)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by