【原创】Matlab教学实验常见问题反馈(by dsj6700417)

dynamic

以下内容由论坛超级版主dsj6700417在Matlab教学实验中，学生问题反馈的总结，在此我们对dsj6700417表示感谢！

我讲一下大家常犯的错误：

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

2.有同学要数值求解多参数方程组，使用了fzero.
  fzero只能求单参数方程的零点.多参数方程的求解用fsolve
    试看下面的例子：

1. x=[0.0350;0.1793;0.5573];%初值
   
2.     options=optimset('MaxFunEvals',30000,'MaxIter',30000);
   
3.     [x,Fval,exitflag] = fzero(@fx,x,options);x
   
4.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%把以下两行百分号之间的部分建立m文件，保存
   
5.     function f=fx(x)
   
6.     F=x(1);K=x(2);K_a=x(3);
   
7.     f=[-62500/3*F*(-exp(-3*K_a)+exp(-3*K))/(-K_a+K)-7639/10;
   
8.      -62500/3*F*(-exp(-18*K_a)+exp(-18*K))/(-K_a+K)-7639/100;
   
9.         -62500/3*F*(-exp(-20*K_a)+exp(-20*K))/(-K_a+K)-267/5];
   
10.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

复制代码

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

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

1.     clear
   
2.     x=0:9;
   
3.     y=[1        2.31969        4.50856        6.90568        6.00512        5.56495        5.32807        7.56101        8.9392        9.5817];
   
4. 
   
5.     f=inline('1/(a(1)+(a(2)-a(1))*exp(-a(3)*x))','a','x');%logstic函数
   
6.     [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时，说明拟合很成功)


我们看一下初值选取要注意的问题，先看下例：

1.     clear
   
2.     t=[ 3   18   20];
   
3.     y=[ 763.9000  76.3900  53.4000];
   
4.     f=inline('a(1)*0.1*(exp(-a(2)*t)-exp(-a(3)*t))/(a(3)-a(2))','a','t');
   
5. 
   
6.     ff=optimset;
   
7.     ff.MaxFunEvals=100000;
   
8.     ff.MaxIter=10000;
   
9.     [xx ,res]=lsqcurvefit(f,[1 ,1 ,1],t,y,[],[],ff)

复制代码

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

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


再看下面的例子(注意 fsolve也属于优化函数)

1. x=[1 1 2];%初值
   
2. [x,Fval,exitflag] = fsolve(@fx,x);x
   
3.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%把以下两行百分号之间的部分建立m文件，保存
   
4.     function f=fx(x)
   
5.     F=x(1);K=x(2);K_a=x(3);
   
6.     f=[-62500/3*F*(-exp(-3*K_a)+exp(-3*K))/(-K_a+K)-7639/10;
   
7.      -62500/3*F*(-exp(-18*K_a)+exp(-18*K))/(-K_a+K)-7639/100;
   
8.         -62500/3*F*(-exp(-20*K_a)+exp(-20*K))/(-K_a+K)-267/5];
   
9.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

复制代码

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

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

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

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

复制代码

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

1. x=[3;4;200];%初值
   
2.     options=optimset('MaxFunEvals',30000,'MaxIter',30000);
   
3.     [x,Fval,exitflag] = fsolve(@fx,x,options);x

复制代码

优化就中止了.


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

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

我把血药浓度问题中，三个参数的求取程序附在后面，提供给大家一个参考，大家有兴趣的可以试着把自己有错误的程序，按照上面的错误警醒，一步步修改一下，体会一下优化需要注意的细节.

1. %%确定F,K,K_a的值
   
2.     x=[0.0350;0.1793;0.5573];%x=[3;4;200];%初值
   
3.     options=optimset('MaxFunEvals',30000,'MaxIter',30000);
   
4.     [x,Fval,exitflag] = fsolve(@fx,x,options);x   % 返回F,K,K_a的具体数值，x=[F,K,K_a]
   
5.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%把以下两行百分号之间的部分建立m文件，保存
   
6.         function f=fx(x)
   
7.         F=x(1);K=x(2);K_a=x(3);
   
8.         f=[-62500/3*F*(-exp(-3*K_a)+exp(-3*K))/(-K_a+K)-7639/10;
   
9.          -62500/3*F*(-exp(-18*K_a)+exp(-18*K))/(-K_a+K)-7639/100;
   
10.             -62500/3*F*(-exp(-20*K_a)+exp(-20*K))/(-K_a+K)-267/5];
    
11.     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

复制代码

lelei0626

这么有用的东西怎么没人支持？   初学者的好东西啊   谢谢楼主  强大的楼主

海上飞鱼

讲的非常好！支持！

seifer08ms

非常实用的内容啊，我支持

付先武

不错,学习了,这么好的东西怎么没人顶啊,我来顶!!!!!!

junjian556

楼主牛人，强烈支持，吼吼

weisinih7osa

这个帖子不回对不起自己！
我觉得MATLAB技术论坛是注册对了！

zhouyang664

支持，支持，再支持！

wgxzy1

强烈支持中..............