lsqcurvefit five parameter estimation with piecewise continuous function

2016 2 16
1 个回答
更新时间：2016 2 19
8 次查看（30 天）

1 个投票

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);

请先登录，再回答此问题。

回答（1 个）

Matt J
Matt J 2016-2-16

0 个投票

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.

5 个评论
显示 3更早的评论 隐藏 3更早的评论

Matthew
Matthew 2016-2-17
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
Matt J 2016-2-17
编辑：Matt J 2016-2-17
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
Matthew 2016-2-19
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
Matthew 2016-2-19
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
Matt J 2016-2-19
Yes, I think it is preferable, but of course you should test it.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Solver Outputs and Iterative Display 的更多信息

标签

提问:

Matthew
Matthew
2016-2-16

评论：

Matt J
Matt J
2016-2-19

Community Treasure Hunt

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

Start Hunting!

Translated by