lsqcurvefit five parameter estimation with piecewise continuous function

I am attempting to use lsqcurvefit to optimize the estimation of five parameters in piecewise continuous function simultaneously, however my script is returning the following error: "Error using symengine. Array sizes must match."

Any comments/suggestions on how to correctly estimate the parameters for this function would be much appreciated.

%%Test parameters and symbolic variables
syms tau
gam=0.01;
t0=11;

%%Simulated stress response to ramp strain and hold experiment
ramp_time=[0,3,9,10.9];
ramp_stress=[0,30,400,455.2];
hold_time=[11.5,23.8,29.8,34.1,38.2,45.2,70.3,84.6,91.3];
hold_stress=[428.6,324.1,267.3,243.2,228.1,213.1,199.6,198.7,198.5];
time=[ramp_time,hold_time];
stress=[ramp_stress,hold_stress];

%%Curve fitting
fun = @(x,time) ((x(1)*x(2)*gam/(1+x(3)*log(x(4)/x(5))))*int((1+x(3)*(igamma(0,(ramp_time-tau)/x(4))-igamma(0,(ramp_time-tau)/x(5))))*exp(x(2)*gam*tau),tau,0,t0)).*(time<t0)+((x(1)*x(2)*gam/(1+x(3)*log(x(4)/x(5))))*int((1+x(3)*(igamma(0,(hold_time-tau)/x(4))-igamma(0,(hold_time-tau)/x(5))))*exp(x(2)*gam*tau),tau,0,t0)).*(time>=t0);
x0=[10,300,300,1,1000];  % Initial guess
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
lb=[];
ub=[];
x = lsqcurvefit(fun,x0,time,stress,lb,ub,options);

Answer 1

[Matt J]

I suggest you forget about lsqcurvefit for the moment and concentrate on getting your objective function to successfully evaluate at x0. It is there that the error is being generated.

[Matthew]

Thank you Matt. That resolved some of the issues, and allowed me to rudimentarily pick a better initial guess. I also reduced the problem by only considering the second half of the objective function, to tackle more manageable parts.

Assuming the objective function is now correct (my initial guess seems realistic and the output of the objective function based off this guess seems realistic), why is the optimization taking so long? When displaying the results of each iteration, I don't seem to successfully evaluate the function even once.

Thank you so much, Matthew.

function QLV_LM

%%Test parameters and symbolic variables
syms tau
gam=0.1;
t0=11;

%%Simulated stress response to ramp strain and hold experiment
hold_time=[11.5,23.8,29.8,34.1,38.2,45.2,70.3,84.6,91.3];
hold_stress=[428.6,324.1,267.3,243.2,228.1,213.1,199.6,198.7,198.5];

%%Curve fitting: Hold Phase Only
x0=[hold_stress(1),1,300,1,1000];  % Initial guess
fun = @(x,hold_time) double((x(1)*x(2)*gam/(1+x(3)*log(x(4)/x(5))))*int((1+x(3)*(vpa(igamma(0,(hold_time-tau)/x(4)))-vpa(igamma(0,(hold_time-tau)/x(5)))))*exp(x(2)*gam*tau),tau,0,t0));

figure(1)
plot(hold_time,hold_stress,'ok',hold_time,fun(x0,hold_time))
xlabel('Time')
ylabel('Stress')
legend('Experiment','Model')

options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt','Display','iter');
lb=[];
ub=[];
[x,resnorm,residual,exitflag,output] = lsqcurvefit(fun,x0,hold_time,hold_stress,lb,ub,options);

[Matt J]

For whatever reason, I'm finding that the time required to evaluate fun(x,hold_time) is very non-repeatable and also can very greatly with x. This may be due to some quirk of the Symbolic Toolbox functions. I don't know if it is wise to use symbolic integration (or symbolic anything else) in conjunction with a numerical minimization routine like lsqcurvefit. If you have the Statistics Toolbox, I imagine that GAMINV might do a better job of those inverse gamma integrals...

[Matthew]

Thank you for the response.

I only partially see where you are going with this... I do have the Statistics Toolbox, however running gaminv(hold_time,0,1) I get a zero, followed by Infs. Is my interpretation wrong?

I know from going through the integration by parts that the breakdown is occurring at the lower bound of integration t=0, which is undefined/INF. Experimentally, I can probably ignore t=0, but that seems like an unjustified assumption on my part regarding the underlying theory.

I am following an example in a text, and the author converts from incomplete gamma functions with upper limits of infinity to a generalized gamma function with upper and lower bounds. I cannot follow his math however, and the text has had errors thus far, so I am not willing to trust it explicitly. It would also be good to thoroughly understand the math/theory. I cannot find any other references to this equation as of yet. For reference, I am attempting to develop a script for Fung's Quasilinear Viscoelastic Theory.

Thank you.

[Matthew]

What do you think about replacing all of the symbolic integrations int() with numeric approximations integral(). Compare the analytical approach

int(exp(-x)/x,x,0,Inf)

to the numerical approach

fun = @(x) exp(-x)./x;
integral(fun,0,Inf)

I don't know why I didn't think of this earlier, but your comment about staying away from symbolic integration made me think of it...

[Matt J]

Yes, I think it is preferable, but of course you should test it.