数值实验三

最优化数值实验 报告 软件系统测试报告下载sgs报告如何下载关于路面塌陷情况报告535n,sgs报告怎么下载竣工报告下载 三 实验目的： 1，​ 能够对具体的问题用合适的最优化方法进行求解。 2，​ 对同一个问题用不同的优化方法进行求解并比较其优劣。 实验内容：（主要是《最优化练习题017》中的3，4，5题） (1)设某实验对象的变化规律可由专业知识得出下列函数 其中a,b,c是反映对象物理特性的待定参数。经过实验,测得数据为 k 1 2 3 4 5 6 7 8 x 0.20 1.0 2.0 3.0 5.0 7 11.0 16.0 y 5.05 8.88 11.63 12.93 14.15 14.73 15.30 15.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=f1^2+f2^2+f3^2+f4^2+f5^2+f6^2+f7^2+f8^2; 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; k x0 min=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+b^2+tan(c)-15; f=f1^2+f2^2+f3^2+f4^2+f5^2+f6^2+f7^2+f8^2+0.5*t*f9^2-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; end x0 min=subs(f1^2+f2^2+f3^2+f4^2+f5^2+f6^2+f7^2+f8^2,[a b c],x0) 结果： 在上述实验中，选择的初始点为（0，0，0），当 ，停止寻优，其中 = 。最后我们取得的最优解为（a,b,c）=(11.4980, -1.0430,4.8220),最优值为0.0284，比上面的情况要差。但是误差可以接受。 我们获得的新的实验对象的变化规律如下： （2） 设某热敏电阻的阻值与温度t和内含某稀土元素钍的含量比q的函数关系为非线性最小二乘问题的数学模型为 其中h,g,m为电阻的物理参数，是模型待定参数。由实验测出下列数据 T 10 20 10 20 1 Q 1.0 1.0 2.0 2.0 0.0 Z 0.126 0.219 0.076 0.126 0.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=f1^2+f2^2+f3^2+f4^2+f5^2; 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; k x0 min=subs(f,x,x0); 结果： 在实验中，选择的初始点为（0，0，0），当 ，停止寻优，其中 = 。最终计算得到的最优解为 。最优值为0.0092。 最终的数学模型可以表示为 这时我们取第二组数据进行验证，发现在 ，和原来的值相差了大概0.06，偏差有些大了。 造成这种现象的原因可能是因为在建立模型的时候出了问题，换成下面的形式就显得合理一些。 这是因为若将约束条件取成上面的除的形式，计算得过程中由于截断误差的原因，会使得结果变得比较差。但若采取这种和差的形式，可能会好些。我们计算得到的数学模型为(其中 = ，因为取得再小的话，在自己电脑上运行时，导致迭代的时间太长。不过虽然精度不是非常高，仍然比上面的情况要好一点。) ， 若写成这种形式，效果要好些。由公式计算得到的值如下 T 10 20 10 20 1 Q 1.0 1.0 2.0 2.0 0.0 Z 0.1026 0.2320 0.0627 0.1349 0.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; k bb0 g1 min=subs(f,b.',bb0); min 结果： 结果感觉并没有完全抑制掉噪声，随机噪声的影响使得每次解出的值都有不小的变化。可能是N取得还不够大的缘故，因为只有当N取得足够大时，matlab产生的随机噪声才比较精确的反应正态随机分布的特性。 实验结论： 通过对上面几个实际问题的研究，觉得在解决实际问题的时候，模型的建立显得非常的重要，就像第二题那种情况一样，要考虑好多因素。模型建好之后，也就是选择合适的方法进行最优化求解了。 附录 进退法确定一维搜索区间 function result=jintuifa(y,r) t0=0; step=0.0125; t1=t0+step; ft0=subs(y,{r},{t0});ft1=subs(y,{r},{t1}); if(ft1<=ft0) step=2*step;t2=t1+step;ft2=subs(y,{r},{t2}); while(ft1>ft2) t1=t2;step=2*step;t2=t1+step;ft1=subs(y,{r},{t1});ft2=subs(y,{r},{t2}); end else step=step/2;t2=t1;t1=t2-step;ft1=subs(y,{r},{t1}); while(ft1>ft0) step=step/2;t2=t1;t1=t2-step;ft1=subs(y,{r},{t1}); end end result=[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)>=e if f1