Estimation of parameters for curvefit using lsqcurvefit

Hi everyone,

On the attached script I am trying to use lsqcurvefit to estimate parameters x(1), x(2), x(3), x(4), x(5) and x(6). This is done in lines 480 to 488, and given as

x0=[8.825e-6;8.4e-3;2e-8;4.7e-3;8.4e-4;9.598e-4];

[x,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(Results,x0,Exp_Time,...
    [Exp_SO2_g;Exp_CO2_g;Exp_S_total;Exp_C_total;Exp_CCaCO3;Exp_CCaSO3;Exp_pH;Exp_pH_Int]);

fprintf(1,'\tParameters:\n')
for k1 = 1:length(x)
    fprintf(1, '\t\tx(%d) = %8.5f\n', k1, x(k1))
end 

I have read/studied the lsqcurve documentation in

https://www.mathworks.com/help/releases/R2018b/optim/ug/lsqcurvefit.html

however, the examples and descriptions given there are too simplistic. I would like to ask for help on how best can i implement lsqcurvefit on my script.

The script run with initial values, however calculated values does not fit the experimental values yet.

Kind Regards

Dursman

Answer 1

[Alan Weiss]

Apparently, you are solving an ODE by forward Euler steps. Don't do that. It is both inaccurate and wasteful of computer time. I refer to these steps in your code:

    dydt = odes (t, y0);
    
    % integrate with explicit Euler
    y = y0 + dydt' * Delt;

Instead, you should call ode45 to solve your ODE.

You also have many more global variables than you need, but that is a discussion for another time.

Now for the question you asked. I get an error when running your supplied code:

Error using lsqfcnchk (line 108)
FUN must be a function, a valid character vector expression, or an inline function object.

Error in lsqnsetup (line 46)
    funfcn = lsqfcnchk(FUN,caller,lengthVarargin,funValCheck,flags.grad);

Error in lsqcurvefit (line 201)
    lsqnsetup(FUN,xCurrent,LB,UB,options,defaultopt,allDefaultOpts, ...

Error in mytry (line 483)
[x,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(Results,x0,Exp_Time,...

Specifically, your Results argument is not a function handle. lsqcurvefit needs to be able to change the parameters and rerun your calculations to fit the function to the data. So I am not completely sure what you are trying to do, but you haven't supplied us with working code to check.

Have you seen the example Fit and Ordinary Differential Equation (ODE)? It is intended to guide you through using lsqcurvefit to fit an ODE.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Dursman Mchabe]

Thanks a lot Alan,

You interpreted my problem accurately. I was forced to use the forward Euler steps by the fact that the 7 ODEs are coupled with 2 algebraic equations, and ODE45 was not suitable. The ODEs are

function dcdt = odes (t, c)
global x V_Headspace F CSO2_in R T HSO2 HCO2 KSPCaSO3  BETCaCO3 MWCaCO3 ...
    KCO2 KHCO3 KSO2 KHSO3 CSO2_bulkslurry  CHSO3 CSO3 CCO2_bulkslurry ...
    CHCO3 CCO3 CCaCO3 CH CCO2_in S_total C_total Ca_total Trial pHInter pg

S_total  = c(3) ;
C_total = c(4) ;
Ca_total = c(5) ;
pH1 = 7; 
pH = fzero (@HionpH, pH1 );
CH = 10^(3-pH) ;
CSO2_bulkslurry     = (S_total*CH^2)/(CH^2 + KSO2*CH + KSO2*KHSO3);
CHSO3               = (S_total*KSO2*CH)/(CH^2 + KSO2*CH + KSO2*KHSO3);
CSO3                = (S_total*KSO2*KHSO3)/(CH^2 + KSO2*CH + KSO2*KHSO3);
CCO2_bulkslurry     = (C_total*CH^2)/(CH^2 + KCO2*CH + KCO2*KHCO3);
CHCO3               = (C_total*KCO2*CH)/(CH^2 + KCO2*CH + KCO2*KHCO3);
CCO3                = (C_total*KCO2*KHCO3)/(CH^2 + KCO2*CH + KCO2*KHCO3);
CH;
pH;
p_exit = c(1)*R*T ;
pg = p_exit;
Trial ;
Tnew = fsolve ( @NSub , Trial);
pHInter = 3-log10(Tnew(3) );
pstar = Tnew(1) * HSO2;
NSO2 = x(4) *(pg-pstar) ;
Trial = Tnew;


%% ODEs
dcdt(1) = (1/V_Headspace) * (F * CSO2_in - F * c(1) ) - NSO2  ;

%-----------------------------------------------------------------
E_CO2 = 1;
NCO2 = x(6)*E_CO2*( c(2)*R*T/HCO2 - CCO2_bulkslurry);
dcdt(2) = (1/V_Headspace) * (F * CCO2_in - F * c(2)) - NCO2 ;
%------------------------------------------------------------------

n= 3;
sat = c(5)*CSO3/ KSPCaSO3  ;
if ( sat < 1.0 )
    RCaSO3  = 0.0; % cut off cacs3 here.
else
    RCaSO3 =   x(3) * ( sat -1)^n ;
end
RCaCO3 = x(1)*BETCaCO3*MWCaCO3*CCaCO3*CH *(1- (x(2)*CH)/(1+x(2)*CH));
dcdt(3) = NSO2 -  RCaSO3;
%---------------------------------------------------------------------

dcdt(4) = NCO2+ RCaCO3;

%---------------------------------------------------------------------

dcdt(5) =  RCaCO3- RCaSO3;

%---------------------------------------------------------------------

if ( c(6) < 1.0e-3 )
    c(6) = 0.0; % cut off caco3 here.
    dcdt(6) = 0.0;
    Time_c = t ; % time when caco3 becomes zero.
    Tag = 1;
else
    dcdt(6) = - RCaCO3 ;
end

%-----------------------------------------------------------------------

dcdt(7) = RCaSO3;

end

Coupled to algebraic equations by

pH1 = 7; 
pH = fzero (@HionpH, pH1 );
CH = 10^(3-pH) ;

and

Trial ;
Tnew = fsolve ( @NSub , Trial);
pHInter = 3-log10(Tnew(3) );

This implies that the output of ODEs (e.g. S_total, Ca_total etc) are used as input for the next timestep in the algebraic equations. The algebraic equations are

function ph =  HionpH (pH)
global KCO2 KHCO3 KSO2 KHSO3 Kw S_total C_total Ca_total

CH =  10^(3-pH);

ph =  CH + 2  * Ca_total - ((S_total*KSO2*CH)/(CH^2 + KSO2*CH + KSO2*KHSO3))...
    - 2*((S_total*KSO2*KHSO3)/(CH^2 + KSO2*CH + KSO2*KHSO3))...
    -((C_total*KCO2*CH)/(CH^2 + KCO2*CH + KCO2*KHCO3))...
    - 2*((C_total*KCO2*KHCO3)/(CH^2 + KCO2*CH + KCO2*KHCO3))- Kw/CH ;
end

and

function ph1 =  HionStar (pH)

global KHCO3 KSO2 KHSO3 Kw Ca_total CHSO3 CSO3 CSO2_inter  CHSO3_inter ...
    CSO3_inter CHCO3_inter  CCO3_inter CHCO3 CCO3

CH =  10^(3-pH);
CSO2_inter ;
CHSO3_inter = KSO2*CSO2_inter  / CH ;
CSO3_inter = KHSO3  *CHSO3_inter  /CH;
Stotal = CHSO3 + CSO3;
Ctotal = CHCO3 + CCO3 ;
CHSO3_inter = Stotal * (   1+  KHSO3  / CH) ;
CSO3_inter = KHSO3  * CHSO3_inter  /CH;
CHCO3_inter = Ctotal * ( 1+  KHCO3  / CH) ;
CCO3_inter = KHCO3  *CHCO3_inter  /CH ;
Ca   = Ca_total;
%---------------------------------------------------------------------------------------
ph1 = CH + 2*Ca -(CHSO3_inter + 2*CSO3_inter + CHCO3_inter + 2* CCO3_inter +  Kw/CH );

end

Similarly, the output of algebraic equations are used as input in the ODEs. The cycle goes on and on for a period of 183 000 seconds.

The code only runs well when the lsqcurvefit steps are commented out.

% x0=[8.825e-6;8.4e-3;2e-8;4.7e-3;8.4e-4;9.598e-4];
% 
% [x,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@SO2_OdeDriver,x0,Exp_Time,...
%     [Exp_SO2_g;Exp_CO2_g;Exp_S_total;Exp_C_total;Exp_CCaCO3;Exp_CCaSO3;Exp_pH;Exp_pH_Int]);
% 
% fprintf(1,'\tParameters:\n')
% for k1 = 1:length(x)
%     fprintf(1, '\t\tx(%d) = %8.5f\n', k1, x(k1))
% end

This is probably because the Results argument is not passed as a function. My wish is to pass the Results argument as a function. I don't know how to do that.

I have looked at Fit and Ordinary Differential Equation (ODE), it does not seem to cater for my current situation.

[Alan Weiss]

I am not 100% sure, but I believe that you can express your problem as a differential-algebraic system of equations. Really, it is best not to use forward Euler steps. Use one of the built-in solvers such as ode15s or ode15i.

To use lsqcurvefit you need to express your results as a function of the input parameters that you are trying to fit. In other words, your 6-D vector x should give you a result that is a (possibly very long) vector of differences between the output of the ODE and the data you are trying to match. In other words, whatever you do, in concept the objective function should look like this:

function resultingdifference = complicatedfunction(x,xdata,otherparameters)

where x is the vector of parameters you are trying to fit, xdata is the data, otherparameters are other problem data that you might want to pass in. The resultingdifference output is a vector predictedresponse - ydata, where ydata is the actual response that you are trying to match.

I hope that this helps orient you.

Alan Weiss

MATLAB mathematical toolbox documentation

[Dursman Mchabe]

Thanks a lot once more Alan. I have been stuck for more than 7 months on this problem. I completely agree that the inbuilt solvers are accurate and reliable. I tried various approaches with the guidence from mathworks members. Ofcourse, I could have been doing it wrongly. But here are some of the thinks I have tried and help/guidence I have received:

https://www.mathworks.com/matlabcentral/answers/403038-is-it-possible-to-optimize-the-parameters-and-initial-conditions-such-that-the-variables-calculated?s_tid=prof_contriblnk

https://www.mathworks.com/matlabcentral/answers/401662-undefined-operator-for-input-arguments-of-type-struct-error-in-decic-line-41-free_yp?s_tid=prof_contriblnk

https://www.mathworks.com/matlabcentral/answers/412248-ode15s-give-an-error-not-enough-input-arguments?s_tid=prof_contriblnk

https://www.mathworks.com/matlabcentral/answers/412270-this-dae-appears-to-be-of-index-greater-than-1?s_tid=prof_contriblnk

https://www.mathworks.com/matlabcentral/answers/414912-error-using-vertcat-dimensions-of-arrays-being-concatenated-are-not-consistent?s_tid=prof_contriblnk

https://www.mathworks.com/matlabcentral/answers/415002-is-there-any-systemic-way-of-determining-a-better-guess-while-solving-daes-using-ode15s?s_tid=prof_contriblnk

https://www.mathworks.com/matlabcentral/answers/415579-plot-coming-up-blank?s_tid=prof_contriblnk

https://www.mathworks.com/matlabcentral/answers/420546-how-can-i-check-which-parameters-is-invalid-in-my-code?s_tid=prof_contriblnk

https://www.mathworks.com/matlabcentral/answers/420847-parameter-estimation-when-solving-dae-using-ode15i-solver?s_tid=prof_contriblnk

https://www.mathworks.com/matlabcentral/answers/421496-estimation-of-parameters-and-initial-conditions-using-lsqcurvefit-gives-an-error-message-index-exc?s_tid=prof_contriblnk

https://www.mathworks.com/matlabcentral/answers/421645-too-many-output-arguments-error-when-using-lsqcurvefit?s_tid=prof_contriblnk

https://www.mathworks.com/matlabcentral/answers/421704-lsqcurvefit-gives-an-error-message-index-exceeds-array-bounds?s_tid=prof_contriblnk

https://www.mathworks.com/matlabcentral/answers/421837-which-input-arguments-are-missing-on-my-code?s_tid=prof_contriblnk

etc.

I am truly stuck, but I'll keep on trying.

[Alan Weiss]

I see that you have been busy. I will give you what advice I can, and hope that it helps in some way.

I think that you need to take a step back and try to figure out what things look like from a high level. Instead of being buried in code, you have a differential-algebraic system that is parametrized by x. You want to find the x that gives a solution as close as possible to some data.

So the highest-level view of your problem is to go from x to the numerical solution of your system, and compare it with the existing data. Once you can have x as an input and the numerical solution as an output, you are ready to feed the problem to lsqcurvefit.

OK, so the next step is to figure out code that goes from x to a solution. You sort of seem to have that as an Euler-stepping monstrosity. Once again, step back and try to formulate the differential equations as a first-order ODE system and the algebraic equations as an associated set of nonlinear equations in the problem variables. Once you have that, you may need to reduce the order of the algebraic system so that the entire system is first order. You can do this. Really! Once you do this, then you won't need Euler steps, but can use a real solver that will converge much, much more accurately and quickly.

So my suggestion, which you might not like, is to first fix your differential-algebraic equations so that they fall into the domain of ode15s and then write your objective function as a short wrapper around the equation solver. Then feed the whole thing to lsqcurvefit.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

[Maryam Alyahyai]

Hi Alan

I have a similar problem. I am not very good in matlab. My problem is about paramater estimation of a system of ODEs. I tried to use the lsqcurvefit in a way similar to fminsearchbnd. I managed to get a good fit and accepted estimation of the parameters when using fminsearchbnd. However, when using lsqcurvefit I ended up with errors!

I have attached the code and the excel sheet for data.

Thanks in advance for your help.

[Alan Weiss]

I suggest that you start a new question and include the entire error messages (everything in red).

Alan Weiss

MATLAB mathematical toolbox documentation

[Maryam Alyahyai]

Done. I've posted it as a new question.

[Maryam Alyahyai]

https://www.mathworks.com/matlabcentral/answers/1693450-fitting-a-system-of-odes-to-data-using-lsqcurvefit-to-estimate-system-s-parameters#comment_2096360

Here is the question.