1、最优化数值实验报告三实验目的：1， 能够对具体的问题用合适的最优化方法进行求解。2， 对同一个问题用不同的优化方法进行求解并比较其优劣。实验内容：（主要是最优化练习题 017中的 3，4，5 题）(1)设某实验对象的变化规律可由专业知识得出下列函数 caexyxb/)(其中a,b,c是反映对象物理特性的待定参数。经过实验,测得数据为k 1 2 3 4 5 6 7 8x 0.20 1.0 2.0 3.0 5.0 7 11.0 16.0y 5.05 8.88 11.63 12.93 14.15 14.73 15.30 15.60请你把此问题表示为最小二乘问题的模型，然后用你认为合适的方法求解，得出
2、待定参数a,b,c的具体数值，将实验对象的变化规律写出。进一步，如果根据专业知识，知道待定参数a,b,c还满足下列条件： ,那么结果又是如何？请你写出求解的数学思15tan2cb想，求解的全过程，并分析你的方法的优缺点。解：首先将问题表示为最小二乘问题的模型，即是将目标函数写成若干个函数的平方和的形式，一般可以写成 mixff12)()(其中 是 中的点。一般假设 ，最小而成问题就是求Tnxx),.(21nRnmixff12)()(i对于本题而言，有m=8,且 60.15)(,30.15)( ,73.14)(,492 68.,.0.16/80.1/7 7/.5.34 0.230./20./1
3、caexfcaexf caexfff bb bbb所以问题转化为求 812)()(minixff通过在前次实验中对各种无约束问题的算法分析和比较知道，BFGS是当中一种性能比较好的，所以首先我们采用BFGS法。（同前两次实验，jintuifa和gold对应进退法和黄金分割法。见附录）程序：syms a b c r;f1=a*exp(b/0.20)+c-5.05;f2=a*exp(b/1.0)+c-8.88;f3=a*exp(b/2.0)+c-11.63;f4=a*exp(b/3.0)+c-12.93;f5=a*exp(b/5.0)+c-14.15;f6=a*exp(b/7.0)+c-14.73
4、;f7=a*exp(b/11.0)+c-15.30;f8=a*exp(b/16.0)+c-15.60;f=f12+f22+f32+f42+f52+f62+f72+f82;x=a,b,c;df=jacobian(f,x);df=df.;error=1e-6;x0=0,0,0;g1=subs(df,x,x0);k=0;H0=1,0 0;0,1 0;0 0 1;while(norm(g1)error)if k=0d=-H0*g1;elseH1=H0+(1+qk*H0*qk/(pk*qk)*(pk*pk)/(pk*qk)-(pk*qk*H0+H0*qk*pk)/(pk*qk);d=-H1*g1;H0=
5、H1;endz=subs(f,x,x0+r*d);result=jintuifa(z,r);result2=gold(z,r,result);step=result2;x0=x0+step*d;g0=g1;g1=subs(df,x,x0);qk=g1-g0; pk=step*d;k=k+1;end;kx0min=subs(f,x,x0)结果：在上述实验中，选择的初始点为（0，0，0），当 ，停止寻优，其中)(1kxf= 。最后我们取得的最优解为（a,b,c）=(11.3457, -1.0730,4.9974),最优值为61*1.9877e-004。效果还是比较不错的。我们可以得到实验对象的变化
6、规律如下： 974.3457.1)(/03.1xexy我们取k=4的时候带入验证， =12.9315，与原来0./.的12.93相差无几。进一步，我们假设参数a,b,c满足下面的条件 ，打算采用拉格朗日15tan2cb乘子法，惩罚函数法和广义乘子法来求解。但是在用拉格朗日乘子法和惩罚函数法求解的时候，发现BFGS法会出现分母为零的情况，导致计算无法正确进行，是由计算过程中的误差的产生而使得分母为0的。但是广义乘子法没有这样的情况。所以最后采用的广义乘子法。程序：syms a b c t v;f1=a*exp(b/0.20)+c-5.05;f2=a*exp(b/1.0)+c-8.88;f3=a*
7、exp(b/2.0)+c-11.63;f4=a*exp(b/3.0)+c-12.93;f5=a*exp(b/5.0)+c-14.15;f6=a*exp(b/7.0)+c-14.73;f7=a*exp(b/11.0)+c-15.30;f8=a*exp(b/16.0)+c-15.60;f9=2*a+b2+tan(c)-15;f=f12+f22+f32+f42+f52+f62+f72+f82+0.5*t*f92-v*f9;x=a,b,c;error1=1e-6;x0=0,0,0;t0=1;c0=1.5;v0=1;j=0;beta=0.5;error2=1e-4;H0=1,0 0;0,1 0;0 0
8、1;k=0;h0=subs(f9,a,b,c,x0(1),x0(2),x0(3);h0=double(h0);while(norm(h0)error1)y=subs(f,t,v,t0,v0);dy=jacobian(y,x);dy=dy.;g1=subs(dy,x,x0);double(g1);while(norm(g1)error2)%BFGS法求解无约束最优化问题if k=0d=-H0*g1;elseH1=H0+(1+qk*H0*qk/(pk*qk)*(pk*pk)/(pk*qk)-(pk*qk*H0+H0*qk*pk)/(pk*qk);d=-H1*g1;H0=H1;endz=subs(y
9、,x,x0+r*d);result=jintuifa(z,r);result2=gold(z,r,result);step=result2;x0=x0+step*d;g0=g1;g1=subs(dy,x,x0);double(g1);qk=g1-g0; pk=step*d;k=k+1;endh1=subs(f9,a,b,c,x0(1),x0(2),x0(3);h1=double(h1);if norm(h1)=beta*norm(h0)t0=t0*c0;endv0=v0-t0*h1;h0=h1;j=j+1;endx0min=subs(f12+f22+f32+f42+f52+f62+f72+f8
10、2,a b c,x0) 结果：在上述实验中，选择的初始点为（0，0，0），当 ，停止寻优，其中)(1kxf= 。最后我们取得的最优解为（a,b,c）=(11.4980, -1.0430,4.8220),最优值为61*0.0284，比上面的情况要差。但是误差可以接受。我们获得的新的实验对象的变化规律如下： 820.44980.1)(/3.1xexy（2） 设某热敏电阻的阻值与温度t和内含某稀土元素钍的含量比q的函数关系为非线性最小二乘问题的数学模型为 mqhtgxtz1),(其中h,g,m为电阻的物理参数，是模型待定参数。由实验测出下列数据T 10 20 10 20 1Q 1.0 1.0 2.0
11、 2.0 0.0Z 0.126 0.219 0.076 0.126 0.186试根据背景建立非线性最小二乘问题的数学模型，并用合适的方法求解，得出模型待定参数，分析你的处理过程，你有什么认识？解：同上题，我们可以建立问题的最小二乘法模型为 412)()(minixff186.0*1)( ,126.0*21)(,7.2 ,9.0,6.0*1)(2 231 mhgxf mhgxfhgxf程序：syms h g m r;f1=h*g*10/(1+h*10+m*1.0)-0.126;f2=h*g*20/(1+h*20+m*1.0)-0.219;f3=h*g*10/(1+h*10+m*2.0)-0.07
12、6;f4=h*g*20/(1+h*20+m*2.0)-0.126;f5=h*g*1/(1+h*1+m*0.0)-0.186;f=f12+f22+f32+f42+f52;x=h,g,m;df=jacobian(f,x);df=df.;error=1e-6;x0=1,1,1;g1=subs(df,x,x0);k=0;H0=1,0 0;0,1 0;0 0 1;while(norm(g1)error)if k=0d=-H0*g1;elseH1=H0+(1+qk*H0*qk/(pk*qk)*(pk*pk)/(pk*qk)-(pk*qk*H0+H0*qk*pk)/(pk*qk);d=-H1*g1;H0=H
13、1;endz=subs(f,x,x0+r*d);result=jintuifa(z,r);result2=gold(z,r,result);step=result2;x0=x0+step*d;g0=g1;g1=subs(df,x,x0);qk=g1-g0; pk=step*d;k=k+1;end;kx0min=subs(f,x,x0);结果：在实验中，选择的初始点为（0，0，0），当 ，停止寻优，其中 = )(1kxf 。最终计算得到的最优解为 。61* )109536.8,062.,28.9),( 717mgh最优值为0.0092。最终的数学模型可以表示为 qtxtz 77109536.81
14、02.94),( 这时我们取第二组数据进行验证，发现在 ，和原来的值相差了大概),(z0.06，偏差有些大了。造成这种现象的原因可能是因为在建立模型的时候出了问题，换成下面的形式就显得合理一些。 412)()(minixff86.0*)1(*)( ,126.0*)2*1(0)(,720 9.2 231hgxf mhgfhf 这是因为若将约束条件取成上面的除的形式，计算得过程中由于截断误差的原因，会使得结果变得比较差。但若采取这种和差的形式，可能会好些。我们计算得到的数学模型为(其中 = ，因为取得再小的话，在自己电脑上运行时，导致迭代的时间太长。210*不过虽然精度不是非常高，仍然比上面的情况
15、要好一点。)，qtxtz4509.273.016),(若写成这种形式，效果要好些。由公式计算得到的值如下T 10 20 10 20 1Q 1.0 1.0 2.0 2.0 0.0Z 0.1026 0.2320 0.0627 0.1349 0.0220第五个参数偏离的比较厉害，其它的偏离偏离大概在0.02左右，但考虑到这时候的精度要比上面小得多 = ，所以有理由相信这种做法是好的。 （精度取的比较大时，210*电脑运行的时间过长。 ）（3） 许多非数学领域的问题最终也常转化为解线性代数方程组的问题。以系统辨识学科为例，它研究如何通过试验测出系统的输入和输出数据，根据这些数据来建立未知系统的数学模型
16、。设某一被辨识系统的输出y(k)与输入u(k)之间的定量关系为（此模型称为MA 模型） )(.)1()(0 nkubkubky其中诸系数为待定的辨识参数。此模型的参数辨识问题是指：已知输出 和输入)(ky诸信息，确定模型系数。)(ku假定在时刻 k=n+1,n+2,L,2n+1 时测量到的数据中不含噪声，那么我们有 )1(.)()(0ubnby221un等一系列方程，系统辨识问题转化为求解线性方程组的问题。可以写成下面的形式。如果输入u(k )不是常量，方程有唯一解。 )12()(,)()2()1()1()()(, 10nybnunuAy n但是，上述处理基于测量数据中没有噪声的假设，显然与实
17、际情况不符，实际上每次测量都包含着难以确定的误差, 每个方程需要加上误差，精确解将不存在了。 yeA为了克服未知噪声的影响，需要增加信息量，比如增加在时刻k=2n+2,2n+3,L,2n+N的测量数据，需要增加信息量，比如增加在时刻k=2n+2,2n+3,L,2n+N 的测量数据，这样我们得到的（1. 6）的系数矩阵A 不是方阵，而是一个“高”的矩阵。出现了一个不相容的矛盾方程组, 这样的线性方程组的解还存在吗？Gauss 提出了最小二乘准则 , 用方程的误差平方和取最小时的 为方程的解。)()(yAeJTT现在我们假设输入为 在取时间间隔 的u值，输出为12sin)(ttut tkt,8/1
18、函数 在上述时刻的离散值。 均为标准随机正态变量。221ln(0)(tty 21,n=10,N=20, 试用最小二乘方法辨识本问题的模型参数，建立本问题的模型。解：问题其实就是一个无约束问题的最优值问题。首先按照题意写出 ,然后对下面的无约yA束问题求最优解。 。Tnb).,(210,)(mi(yA其中tktkt tyeku NnyNnuNnunA |)1l(0)(,|)si() )2()1(,12(2)2()(1 221 30(30其中 可由函数randn生成。21,模型建立完毕。（下面给出代码，代码中描述了整个思想。）程序：syms t b0 b1 b2 b3 b4 b5 b6 b7 b8
19、 b9 b10 r;ut=exp(-2*t)*sin(pi*t);yt=10*log(t*t-1);n=10;N=20;delt=1/8;u=zeros(2*n+N,1);y=zeros(n+N,1);A=zeros(n+N,n+1);for i=1:2*n+N %得到u矩阵u(i)=subs(ut,t,i*delt);endfor j=(n+1):(2*n+N) %得到y矩阵y(j-n)=subs(yt,t,j*delt);endu=u+randn(2*n+N,1);y=y+randn(n+N,1);m=n+1; %叠加上高斯噪声，计算A矩阵for i=1:N+nfor j=1:(n+1)A
20、(i,j)=u(m+1-j);endm=m+1;endb=b0 b1 b2 b3 b4 b5 b6 b7 b8 b9 b10.;%最优值求解f=(A*b-y).*(A*b-y);df=jacobian(f,b);df=df.;k=0;error=1e-6;bb0=zeros(n+1,1);g1=subs(df,b.,bb0);k=0;H0=eye(n+1);while(norm(g1)error)if k=0d=-H0*g1;elseH1=H0+(1+qk*H0*qk/(pk*qk)*(pk*pk)/(pk*qk)-(pk*qk*H0+H0*qk*pk)/(pk*qk);d=-H1*g1;H0
21、=H1;endz=subs(f,b.,bb0+r*d);result=jintuifa(z,r);result2=gold(z,r,result);step=result2;bb0=bb0+step*d;g0=g1;g1=subs(df,b.,bb0);qk=g1-g0; pk=step*d;k=k+1;end;kbb0g1min=subs(f,b.,bb0);min结果：结果感觉并没有完全抑制掉噪声，随机噪声的影响使得每次解出的值都有不小的变化。可能是N取得还不够大的缘故，因为只有当N取得足够大时，matlab产生的随机噪声才比较精确的反应正态随机分布的特性。实验结论：通过对上面几个实际问题
22、的研究，觉得在解决实际问题的时候，模型的建立显得非常的重要，就像第二题那种情况一样，要考虑好多因素。模型建好之后，也就是选择合适的方法进行最优化求解了。附录进退法确定一维搜索区间function result=jintuifa(y,r)t0=0; step=0.0125;t1=t0+step;ft0=subs(y,r,t0);ft1=subs(y,r,t1);if(ft1ft2)t1=t2;step=2*step;t2=t1+step;ft1=subs(y,r,t1);ft2=subs(y,r,t2);endelse step=step/2;t2=t1;t1=t2-step;ft1=subs(
23、y,r,t1);while(ft1ft0)step=step/2;t2=t1;t1=t2-step;ft1=subs(y,r,t1);endendresult=t2;黄金分割法进行一维搜索function result=gold(y,r,m)a=0;b=m;e=1e-5;a1=a+0.382*(b-a);f1=subs(y,r,a1);a2=a+0.618*(b-a);f2=subs(y,r,a2);while abs(b-a)=eif f1f2b=a2;a2=a1;f2=f1;a1=a+0.382*(b-a);f1=subs(y,r,a1);elsea=a1;a1=a2;f1=f2;a2=a+0.618*(b-a);f2=subs(y,r,a2);endend;answer=0.5*(a+b);result=answer;
