Fitting model to titration data

I'm trying to fit a model of titration to experimental data. The function "calcpH" calculate a pH value for each experimental point and I'm trying to obtain the value of Ca (and eventually Ka) that fit best the model to the experimental values.

When I uncomment the lines of the fitting procedure I get the following error:

Error using ajuste_listado_1>sseval
Too many output arguments.
Error in ajuste_listado_1>@(x)sseval(x,vbr,pHr)
Error in fminsearch (line 200)
fv(:,1) = funfcn(x,varargin{:});

What is going wrong ????

THANK YOU

CODE OF THE SCRIPT

clear
clc
pHr=[3.15,3.30,3.48,3.65,3.73,3.86,3.93,4.03,4.10,4.20,4.28,4.35,4.44,4.51,4.59,4.66,4.75,4.84,4.94,5.03,5.18,5.34,5.54,5.91,6.10,6.40,8.66,9.32,9.88,10.16,10.39,10.49,10.64,10.97,11.16,11.28,11.37,11.45,11.51,11.56,11.60,11.67];
vbr=[0.00 ,0.20 ,0.40 ,0.60 ,0.80 ,1.00 ,1.20 ,1.40 ,1.60 ,1.80 ,2.00 ,2.20 ,2.40 ,2.60 ,2.80 ,3.00 ,3.20 ,3.40 ,3.60 ,3.80 ,4.00 ,4.20 ,4.40 ,4.60 ,4.70 ,4.80 ,4.90 ,5.00 ,5.10 ,5.20 ,5.30 ,5.40 ,5.50 ,6.00 ,6.50 ,7.00 ,7.50 ,8.00 ,8.50 ,9.00 ,9.50 ,10.00];
Ka=1.8e-5;
Ca=0.01
Va=50;
Cb=0.1;
Vb=0:0.2:10;
Kw=1e-14;

% FITTING 
% fun=@(x)sseval(x,vbr,pHr);
% bestx = fminsearch(fun,x0)
% Ca = bestx(1);

x=zeros(1)
pH = calcpH(x,Ka,Ca,Kw,Cb,Va,vbr);

plot(vbr,pHr,'o',vbr,pH,'x')
title('Titration of acetic acid')
xlabel('Vol NaOH [mL]')
ylabel('pH')

function sseval(x,vbr,pHr)
Ca=x(1)
pH = calcpH(x,Ka,Ca,Kw,Cb,Va,vbr)
sse = sum((pHr - pH').^2)
end

function pH=calcpH(x,Ka,Ca,Kw,Cb,Va,vbr)
VT=Va+vbr;
for i=1:length(vbr)
syms xh
h=vpasolve(Cb*vbr(i)/VT(i)-(Ca*Va/VT(i))*Ka/(Ka+xh)-Kw/xh+xh==0,xh);
pH(i)= -log10(h(3));
end
end

Answer 1

[Rafael Ibáñez]

Thank you very much

Answer 2

[David Wilson]

For one thing, you don't assign an answer to function sseval. That's not going to help things.

Here's my solution:

I've re-written your internal function and removed the numerical root finder. It's only a cubic, so we should just use roots. HOWEVER you have assumed you only get one physcially sensible (i.e. positive) root and use that. Are you sure that's OK? (I followed your approach, so I'm counting on that.)

function pH=calcpH(Ca,Ka,Kw,Cb,Va,vbr)

VT=Va+vbr;

a = Cb*vbr./VT;
b = Ca*Va*Ka./VT; 
c = Ka; 
d = Kw; 
pH = NaN(size(vbr)); % placeholder 

for i=1:length(vbr)
   coeff3 = [1, (a(i)+c), -(b(i)+d-a(i)*c), -c*d]; 
   h = sort(roots(coeff3)); % assume last is the only (sensible) positive  
   h = h(end); 
   pH(i,1)= -log10(h);
end 

end

Now I wrap the objective function around that. I'm using the 2-norm, nice and well-conditioned. The output of this function is ideally small, or better still zero.

function sse = sseval(Ca,vbr,pHr, ...
                      Ka,Kw,Cb,Va)

pH=calcpH(Ca,Ka,Kw,Cb,Va,vbr); 
sse = norm(pHr-pH);  % sum((pHr - pH).^2);
end

Note that you need to pass through all your constants etc correctly.

Now we are ready to call the optimiser from the top script. I like to make sure everything is in columns (your original version had a mixture which was just asking for trouble.) You need a guess for Ca.

pHr = pHr(:); vbr = vbr(:); % ensure columns 

fun=@(x) sseval(x, vbr,pHr, Ka,Kw,Cb,Va);
Ca_guess = 0.01; 

optns = optimset('fminsearch');
[Ca, fval, exitflag] = fminsearch(fun,Ca_guess, optns)

%Ca = fminunc(fun,Ca_guess)

pH = calcpH(Ca,Ka,Kw,Cb,Va,vbr);

plot(vbr,pHr,'o',vbr,pH,'-')
title('Titration of acetic acid')
xlabel('Vol NaOH [mL]')
ylabel('pH')

Finally we should plot the comparison. HOWEVER mine's not that wonderful. I'd hoped for better.

If you want a better fit, you need to "adjust" your constants. For example, Kw = 1.8e-14 gives a much better fit. (Note Kw is a function of temperature, and 1e-14 is not a golden constant! Similarly for your other constants.

[Rafael Ibáñez]

Thnaks¡¡¡¡

It works