数学建模社区-数学中国

标题: Matlab教学实验常见问题反馈(by dsj6700417) [打印本页]

作者: 建不了的模。 时间: 2014-9-18 14:24
标题: Matlab教学实验常见问题反馈(by dsj6700417)
我讲一下大家常犯的错误：

1.有同学要拟合的目标函数是a*x+b，使用了最小二乘拟合lsqcurvefit，或者非线拟合nlinfit函数，
  应该使用多项式拟合polyfit.具体到此例使用polyfit(x,y,1)

2.有同学要数值求解多参数方程组，使用了fzero.
  fzero只能求单参数方程的零点.多参数方程的求解用fsolve
试看下面的例子：
x=[0.0350;0.1793;0.5573];%初值
options=optimset('MaxFunEvals',30000,'MaxIter',30000);
[x,Fval,exitflag] = fzero(@fx,x,options);x
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%把以下两行百分号之间的部分建立m文件，保存
function f=fx(x)
F=x(1);K=x(2);K_a=x(3);
f=[-62500/3*F*(-exp(-3*K_a)+exp(-3*K))/(-K_a+K)-7639/10;
   -62500/3*F*(-exp(-18*K_a)+exp(-18*K))/(-K_a+K)-7639/100;
      -62500/3*F*(-exp(-20*K_a)+exp(-20*K))/(-K_a+K)-267/5];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
复制代码

上面的程序格式上是正确的，但报错的原因在于错误使用了fzero，改为fsolve即可运行

3.我以lsqcurvefit为例，讲一下这个函数使用上的技巧，并解释大家的疑问
先讲使用lsqcurvefit的四个注意事项：
(1)注意模型函数的正确性，
(2)注意初值的选择
(3)注意optimset的设置
(4)注意res的大小，和exitflag的返回值

  我先用下面的代码分析一下：
clear
x=0:9;
y=[1       2.31969       4.50856       6.90568       6.00512       5.56495       5.32807       7.56101       8.9392       9.5817];

f=inline('1/(a(1)+(a(2)-a(1))*exp(-a(3)*x))','a','x');%logstic函数
[xx ,res]=lsqcurvefit(f,[1 1 1],x,y)
复制代码
这是lsqcurvefit最常用的使用方法.上面程序的报错原因在于，inline函数里面的分母是个数组，数组除法要用点除，
在1后面加个点，即1./(a(1)+(a(2)-a(1))*exp(-a(3)*x))，程序就可以运行了.

修改后的程序，残差res很大，任凭你怎么设置optimset的精度及容差，或者修改初值，这个残差res都很大，这时需要考虑函数模型的正确性。logstic函数是个单调函数，而数据的散点图不是单调的，中间一段有下降部分，这个是导致res过大的原因，所以兔子问题我们根据logstic函数的单调性分为三段拟合。

我把问题提的尖锐一些：你不可能用指数函数exp(k^2*x),去拟合散点图呈正弦分布的数据，而要求残差很小；你不可能用一个恒正的函数模型去拟合全负的数据;你又不能用一个单调增函数去拟合一组单调减的数据。

后面这种情况，兔子问题的第二段曲线拟合就是个例子，初值都取正数时，logstic函数就是单调增函数，这样去拟合一组单调减的数据，结果是很差的，解决办法是改最后一个初值为负.如果你细心分析的话，会发现第三个初值为负就可以保证logstic函数是单调减的。

这个例子说明：
(1)要保证函数模型正确，不能使用错误的或者不合适的函数模型；
(2)构建函数时，要注意数组运算符的正确使用；
(3)初值的设置与函数模型的性质有关，要学会大体分析函数模型的性质.拟合残差过大时，试着多改几次初值；
(4)注意观察残差res的大小，res越小，拟合效果越好；
(如果使用[x,resnorm,residual,exitflag]=lsqcurvefit(f,[1 1 1],x,y) ，则exitflag返回-1时，说明拟合很成功)

我们看一下初值选取要注意的问题，先看下例：
clear
t=[ 3 18 20];
y=[ 763.9000  76.3900  53.4000];
f=inline('a(1)*0.1*(exp(-a(2)*t)-exp(-a(3)*t))/(a(3)-a(2))','a','t');

ff=optimset;
ff.MaxFunEvals=100000;
ff.MaxIter=10000;
[xx ,res]=lsqcurvefit(f,[1 ,1 ,1],t,y,[],[],ff)
复制代码

程序报错的原因在于inline函数的最后：(a(3)-a(2))在取初值时等于了0，导致分母为零，优化直接中止.一般情况下，大家习惯取初值为全1矩阵.这也要视具体函数而定.

我们把初值改为[.1,.1,.2]，程序即可运行，但此后不管我们怎么修改初值(要保证后两个不相等)，残差res都大的惊人，这时的原因就在于函数模型给的不对(不是指函数格式不正确，是指问题的函数模型有误)，或者说是函数模型使用的数据不对，上例正是这个原因：函数模型错误

再看下面的例子(注意 fsolve也属于优化函数)
x=[1 1 2];%初值
[x,Fval,exitflag] = fsolve(@fx,x);x
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%把以下两行百分号之间的部分建立m文件，保存
function f=fx(x)
F=x(1);K=x(2);K_a=x(3);
f=[-62500/3*F*(-exp(-3*K_a)+exp(-3*K))/(-K_a+K)-7639/10;
   -62500/3*F*(-exp(-18*K_a)+exp(-18*K))/(-K_a+K)-7639/100;
      -62500/3*F*(-exp(-20*K_a)+exp(-20*K))/(-K_a+K)-267/5];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
复制代码

上面初值的后面两个不相等也是为了避免函数分母为零.

为了得到优化结果,好多人都会用下面的办法：

随意设置一个初值，把运行结果作为初值，反复执行，直到初值和结果相等，即迭代至不动点为止。这个办法是很耗功耗时的，有一个一劳永逸的办法是：在optimset里面调大函数调用次数MaxFunEvals和迭代次数MaxIter，但很多情况上面的办法都不起效果，比如上面的程序改为：

  x=[1 1 2];%初值
options=optimset('MaxFunEvals',30000,'MaxIter',30000);
[x,Fval,exitflag] = fsolve(@fx,x,options);x
复制代码

不管你把MaxFunEvals和MaxIter的值调的再大，结果都不会收敛。这时必须更换初值的类型，比如一个好的运气，把初值改为

x=[3;4;200];%初值
options=optimset('MaxFunEvals',30000,'MaxIter',30000);
[x,Fval,exitflag] = fsolve(@fx,x,options);x
复制代码
优化就中止了.

总结一下，优化时综合考虑以下问题：
(1)函数模型的正确性，函数的正确书写
(2)初值的选择，主要与模型函数的性质有关
(3)optimset的设置，主要用于设置精度，容差，最大调用次数，最大迭代次数
(4)学会查看优化效果：res和exitflag

优化问题本来就是很复杂的，要想得到需要的结果需要设置很多选项和考虑很多因素，鉴于此，不提倡大家使用cftool!

我把血药浓度问题中，三个参数的求取程序附在后面，提供给大家一个参考，大家有兴趣的可以试着把自己有错误的程序，按照上面的错误警醒，一步步修改一下，体会一下优化需要注意的细节.
%%确定F,K,K_a的值
x=[0.0350;0.1793;0.5573];%x=[3;4;200];%初值
options=optimset('MaxFunEvals',30000,'MaxIter',30000);
[x,Fval,exitflag] = fsolve(@fx,x,options);x % 返回F,K,K_a的具体数值，x=[F,K,K_a]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%把以下两行百分号之间的部分建立m文件，保存
      function f=fx(x)
      F=x(1);K=x(2);K_a=x(3);
      f=[-62500/3*F*(-exp(-3*K_a)+exp(-3*K))/(-K_a+K)-7639/10;
      -62500/3*F*(-exp(-18*K_a)+exp(-18*K))/(-K_a+K)-7639/100;
         -62500/3*F*(-exp(-20*K_a)+exp(-20*K))/(-K_a+K)-267/5];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
复制代码

作者: 深V礼 时间: 2014-9-18 17:26
顶顶，好样的

欢迎光临数学建模社区-数学中国 (http://www.madio.net/)