最优化数值实验三.doc

最优化数值实验报告三 实验目的：1， 能够对具体的问题用合适的最优化方法进行求解。2， 对同一个问题用不同的优化方法进行求解并比较其优劣。实验内容：（主要是最优化练习题017中的3，4，5题）(1)设某实验对象的变化规律可由专业知识得出下列函数其中a,b,c是反映对象物理特性的待定参数。经过实验,测得数据为k12345678x0.201.02.03.05.0711.016.0y5.058.8811.6312.9314.1514.7315.3015.60请你把此问题表示为最小二乘问题的模型，然后用你认为合适的方法求解，得出待定参数a,b,c的具体数值，将实验对象的变化规律写出。进一步，如果根据专业知识，知道待定参数a,b,c还满足下列条件：,那么结果又是如何？请你写出求解的数学思想，求解的全过程，并分析你的方法的优缺点。解：首先将问题表示为最小二乘问题的模型，即是将目标函数写成若干个函数的平方和的形式，一般可以写成其中是中的点。一般假设，最小而成问题就是求对于本题而言，有m=8,且所以问题转化为求通过在前次实验中对各种无约束问题的算法分析和比较知道，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;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=0 d=-H0*g1; else H1=H0+(1+qk*H0*qk/(pk*qk)*(pk*pk)/(pk*qk)-(pk*qk*H0+H0*qk*pk)/(pk*qk); d=-H1*g1; H0=H1; end z=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），当，停止寻优，其中= 。最后我们取得的最优解为（a,b,c）=(11.3457, -1.0730,4.9974),最优值为1.9877e-004。效果还是比较不错的。 我们可以得到实验对象的变化规律如下： 我们取k=4的时候带入验证，=12.9315，与原来的12.93相差无几。进一步，我们假设参数a,b,c满足下面的条件，打算采用拉格朗日乘子法，惩罚函数法和广义乘子法来求解。但是在用拉格朗日乘子法和惩罚函数法求解的时候，发现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*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 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=0 d=-H0*g1; else H1=H0+(1+qk*H0*qk/(pk*qk)*(pk*pk)/(pk*qk)-(pk*qk*H0+H0*qk*pk)/(pk*qk); d=-H1*g1; H0=H1; end z=subs(y,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; end h1=subs(f9,a,b,c,x0(1),x0(2),x0(3);h1=double(h1); if norm(h1)=beta*norm(h0) t0=t0*c0; end v0=v0-t0*h1; h0=h1; j=j+1;endx0min=subs(f12+f22+f32+f42+f52+f62+f72+f82,a b c,x0) 结果：在上述实验中，选择的初始点为（0，0，0），当，停止寻优，其中= 。最后我们取得的最优解为（a,b,c）=(11.4980, -1.0430,4.8220),最优值为0.0284，比上面的情况要差。但是误差可以接受。我们获得的新的实验对象的变化规律如下：（2） 设某热敏电阻的阻值与温度t和内含某稀土元素钍的含量比q的函数关系为非线性最小二乘问题的数学模型为其中h,g,m为电阻的物理参数，是模型待定参数。由实验测出下列数据T102010201Q1.01.02.02.00.0Z0.1260.2190.0760.1260.186试根据背景建立非线性最小二乘问题的数学模型，并用合适的方法求解，得出模型待定参数，分析你的处理过程，你有什么认识？解：同上题，我们可以建立问题的最小二乘法模型为程序：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.076;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=0 d=-H0*g1; else H1=H0+(1+qk*H0*qk/(pk*qk)*(pk*pk)/(pk*qk)-(pk*qk*H0+H0*qk*pk)/(pk*qk); d=-H1*g1; H0=H1; end z=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），当，停止寻优，其中= 。最终计算得到的最优解为。最优值为0.0092。最终的数学模型可以表示为这时我们取第二组数据进行验证，发现在，和原来的值相差了大概0.06，偏差有些大了。造成这种现象的原因可能是因为在建立模型的时候出了问题，换成下面的形式就显得合理一些。 这是因为若将约束条件取成上面的除的形式，计算得过程中由于截断误差的原因，会使得结果变得比较差。但若采取这种和差的形式，可能会好些。我们计算得到的数学模型为(其中= ，因为取得再小的话，在自己电脑上运行时，导致迭代的时间太长。不过虽然精度不是非常高，仍然比上面的情况要好一点。)，若写成这种形式，效果要好些。由公式计算得到的值如下T102010201Q1.01.02.02.00.0Z0.10260.23200.06270.13490.0220第五个参数偏离的比较厉害，其它的偏离偏离大概在0.02左右，但考虑到这时候的精度要比上面小得多= ，所以有理由相信这种做法是好的。（精度取的比较大时，电脑运行的时间过长。）（3） 许多非数学领域的问题最终也常转化为解线性代数方程组的问题。以系统辨识学科为例，它研究如何通过试验测出系统的输入和输出数据，根据这些数据来建立未知系统的数学模型。设某一被辨识系统的输出y(k)与输入u(k)之间的定量关系为（此模型称为MA 模型）其中诸系数为待定的辨识参数。此模型的参数辨识问题是指：已知输出和输入诸信息，确定模型系数。假定在时刻 k=n+1,n+2,L,2n+1 时测量到的数据中不含噪声，那么我们有等一系列方程，系统辨识问题转化为求解线性方程组的问题。可以写成下面的形式。如果输入u(k)不是常量，方程有唯一解。但是，上述处理基于测量数据中没有噪声的假设，显然与实际情况不符，实际上每次测量都包含着难以确定的误差, 每个方程需要加上误差，精确解将不存在了。为了克服未知噪声的影响，需要增加信息量，比如增加在时刻k=2n+2,2n+3,L,2n+N的测量数据，需要增加信息量，比如增加在时刻k=2n+2,2n+3,L,2n+N 的测量数据，这样我们得到的（1. 6）的系数矩阵A 不是方阵，而是一个“高”的矩阵。出现了一个不相容的矛盾方程组, 这样的线性方程组的解还存在吗？Gauss 提出了最小二乘准则, 用方程的误差平方和取最小时的为方程的解。现在我们假设输入为在取时间间隔的u值，输出为函数在上述时刻的离散值。均为标准随机正态变量。n=10,N=20, 试用最小二乘方法辨识本问题的模型参数，建立本问题的模型。解： 问题其实就是一个无约束问题的最优值问题。首先按照题意写出,然后对下面的无约束问题求最优解。,其中 其中可由函数randn生成。模型建立完毕。（下面给出代码，代码中描述了整个思想。）程序：syms t b0 b1 b2 b3 b4 b5 b6 b7 b8 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);end for j=(n+1):(2*n+N) %得到y矩阵 y(j-n)=subs(yt,t,j*delt);end u=u+randn(2*n+N,1);y=y+randn(n+N,1);m=n+1; %叠加上高斯噪声，计算A矩阵for i=1:N+n for j=1:(n+1) A(i,j)=u(m+1-j); end m=m+1;end b=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=0 d=-H0*g1; else H1=H0+(1+qk*H0*qk/(pk*qk)*(pk*pk)/(pk*qk)-(pk*qk*H0+H0*qk*pk)/(pk*qk); d=-H1*g1; H0=H1; end z=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产生的随机噪声才比较精确的反应正态随机分布的特性。实验结论： 通过对上面几个实际问题的研究，觉得在解决实际问题的时候，模型的建立显得非常的重要，就像第二题那种情况一样，要考虑好多因素。