24小时热门版块排行榜

返回列表

anarky

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 150
帖子: 16
在线: 8.2小时
虫号: 2655893
注册: 2013-09-15
专业: 机械结构强度学

[求助] 请教matlab非线性拟合结果不唯一的问题已有2人参与

我用matlab拟合一组实验数据，函数是y=a*（b+x)^c+d*(1-f*exp(-g*x)),现在的问题是发现只要改一次初值就会得到不同的拟合结果，或者改一下上下限，也会得到不同的拟合结果，这么一来，我觉得就不可能找到最优解，虽然上下限可以取得很大，认为是无界，但是初值怎么取？因为我不可能尝试所有的参数初值，因为我要拟合多组实验数据，获得公式中参数a(1)~a(6)的变化规律，如果参数结果不唯一，我怎么确定哪个参数才是最优的参数呢？事实上，如果从物理意义来限定参数取值也可以，但是那样得到的参数并不是最优解，而最优解的参数又不一定符合物理意义。我该怎么处理呢？哪位能帮忙看一下，非常感谢！
实验数据和代码如下：
主程序
clc
xdata=load('78028-ta.txt');  %read section1
size1=size(xdata);             %section1 行、列数
sizen(1)=size1(1);             %section1数据量
for i=1:size1;
xd(i)=xdata(i,1);
yd(i)=xdata(i,2);
end

set(gca,'Fontsize',15);
plot(xd,yd,'ro');
xlabel('True Strain');
ylabel('True Stress/MPa');

hold on
lb=[-100000, 0, 0, -10000, -10000, -10000];
ub=[ 100000, 0.02, 1,  10000,  10000,  10000];
a1=1898;a2=0.00002;a3=0.136;a4=542;a5=0.784;a6=62.4;
a0=[a1,a2,a3,a4,a5,a6];
options = optimset('TolFun',1e-6,'LargeScale','on','TolCon',1.0e-6,'MaxFunEvals',2.e10,'MaxIter',2000000,'TolX',1e-6);
exitflag=0;
while (exitflag<=0)                               %反复迭代
[a,resnorm,residual,exitflag]=lsqnonlin(@SV,a0,lb,ub,options);    %nonlinear fitting, solve coefficient a
a0=a;
end
R2=1-resnorm/sum((yd-mean(yd)).^2);
standard=sqrt(resnorm/size(residual,2));
stp=1;
out(stp,1)=R2;
out(stp,2)=standard;
out(stp,3)=a(1);
out(stp,4)=a(2);
out(stp,5)=a(3);
out(stp,6)=a(4);
out(stp,7)=a(5);
out(stp,8)=a(6);

hold on
fid=fopen('SV-fit-78028-ta.txt','wt');
for i=1:size1;
   G(i)=a(1)*(a(2)+xd(i))^a(3)+a(4)*(1-a(5)*exp(-a(6)*xd(i)));
fprintf(fid,'%8.6f %10.6f\n',xd(i),G(i));
end

fprintf(fid,'    R2       standard    s1          s2          s3       V1       V2       V3\n');
fprintf(fid,'%10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f\n',out(i,1),out(i,2),out(i,3),out(i,4),out(i,5),out(i,6),out(i,7),out(i,8));
fclose(fid);
plot(xd,G,'b','LineWidth',2)
hold on
legend('EXP','SV');
hold off
clear
函数
function F=SV(a)
xdata=load('78028-ta.txt');
size1=size(xdata);
sizen(1)=size1(1);

for i=1:sizen(1);
xd(i)=xdata(i,1);
yd(i)=xdata(i,2);
end

for i=1:sizen(1);
F(i)=a(1)*(a(2)+xd(i))^a(3)+a(4)*(1-a(5)*exp(-a(6)*xd(i)));
F(i)=F(i)-yd(i);
end

实验数据
78028-ta.txt
0 223.634
0.0005 261.9
0.001 373.358
0.0015 425.658
0.002 462.61
0.0025 492.699
0.003 513.549
0.0035 524.142
0.004 538.681
0.0045 546.22
0.005 554.575
0.0055 563.233437
0.006 570.27162
0.0065 577.07317
0.007 583.656651
0.0075 590.037892
0.008 596.230537
0.0085 602.246461
0.009 608.09609
0.0095 613.788651
0.01 619.332371
0.0105 624.734633
0.011 630.002102
0.0115 635.14083
0.012 640.156342
0.0125 645.053702
0.013 649.837576
0.0135 654.512279
0.014 659.081814
0.0145 663.549911
0.015 667.920058
0.0155 672.19552
0.016 676.379368
0.0165 680.474496
0.017 684.483635
0.0175 688.409373
0.018 692.254165
0.0185 696.020341
0.019 699.710123
0.0195 703.325629
0.02 706.4514911
0.025 735.7217232
0.03 758.4120509
0.035 778.3482998
0.04 794.7567166
0.045 809.6868333
0.05 822.8422861
0.055 834.8128986
0.06 846.1687589
0.065 856.1567004
0.07 866.0042315
0.075 874.4658077
0.08 882.1951827
0.085 890.0445724
0.09 897.0900675
0.095 903.7413962
0.1 910.0615411
0.105 915.8285085
0.11 921.8803708
0.115 927.2834332
0.12 932.4186891
0.125 937.3558307
0.315 1073.901289

回复此楼

1楼 2015-12-07 21:45:34

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

dingd

铁杆木虫 (职业作家)

计算强帖: 4
应助: 1641 (讲师)
金币: 15037.3
散金: 101
红花: 234
帖子: 3410
在线: 1223.7小时
虫号: 291104
注册: 2006-10-28

【答案】应助回帖

★ ★ ★ ★
感谢参与，应助指数 +1
anarky: 金币+4, ★★★很有帮助, 非常感谢！ 2015-12-09 21:24:24

用1stOpt，不用初值，结果唯一稳定：

CODE:

Function  y=a*(b+x)^c+d*(1-f*exp(-g*x));

Data;

0        223.634

0.0005        261.9

0.001        373.358

0.0015        425.658

0.002        462.61

0.0025        492.699

。。。。。

均方差(RMSE):6.83618029811499
残差平方和(SSE):2944.20174730514
相关系数(R): 0.99915354759785
相关系数之平方(R^2): 0.998307811677368
确定系数(DC): 0.998307811677368
卡方系数(Chi-Square): 4.49199298531208
F统计(F-Statistic): 6736.83883025038

参数最佳估算
---------- -------------
a -1838.97884865162
b 0.00549570600631707
c -0.0663243340918725
d 3047.85332211066
f 0.0810032510811905
g 817.309082573211

其它软件很难得到上面的结果