最优化数值实验报告三

最优化数值实验 报告 软件系统测试报告下载sgs报告如何下载关于路面塌陷情况报告535n,sgs报告怎么下载竣工报告下载 三 实验目的: 1， 能够对具体的问题用合适的最优化方法进行求解。 2， 对同一个问题用不同的优化方法进行求解并比较其优劣。 实验内容:(主要是《最优化 练习题 用券下载整式乘法计算练习题幼小衔接专项练习题下载拼音练习题下载凑十法练习题下载幼升小练习题下载免费 017》中的3，4，5题) (1)设某实验对象的变化规律可由专业知识得出下列函数 b/xy(x),ae，c 其中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的具体数值，将实验对象的变化规律写出。进一步，如果根据专业知识，知道待定参 2数a,b,c还满足下列条件:,那么结果又是如何,请你写出求解的数学思2a，b，tanc,15 想，求解的全过程，并分析你的方法的优缺点。 解: 首先将问题表示为最小二乘问题的模型，即是将目标函数写成若干个函数的平方和的形式，一般可以写成 m2f(x),f(x) ,i,1i nTm,n其中是中的点。一般假设，最小而成问题就是求 Rx,(x,x,...,x)12n m2minf(x),f(x) ,i,1i 对于本题而言，有m=8,且 b/0.20b/1.0b/2.0f(x),ae，c,5.05,f(x),ae，c,8.88,f(x),ae，c,11.63,123 b/3.0b/5.0b/7f(x),ae，c,12.93,f(x),ae，c,14.15,f(x),ae，c,14.73, 456 b/11.0b/16.0f(x),ae，c,15.30,f(x),ae，c,15.6078 所以问题转化为求 82minf(x),f(x) ,ii,1 通过在前次实验中对各种无约束问题的算法分析和比较知道，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) 结果: k，1,,f(x),, 在上述实验中，选择的初始点为(0，0，0)，当，停止寻优，其中= ,6。最后我们取得的最优解为(a,b,c)=(11.3457, -1.0730,4.9974),最优值为1*10 1.9877e-004。效果还是比较不错的。 我们可以得到实验对象的变化规律如下: ,1.0730/xy(x),11.3457e，4.9974 ,1.0730/3.0y(3.0),11.3457e，4.9974 我们取k=4的时候带入验证，=12.9315，与原来的 12.93相差无几。 2进一步，我们假设参数a,b,c满足下面的条件，打算采用拉格朗日乘2a，b，tanc,15 子法，惩罚函数法和广义乘子法来求解。但是在用拉格朗日乘子法和惩罚函数法求解的时候， 发现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) 结果: k，1,,f(x),,在上述实验中，选择的初始点为(0，0，0)，当，停止寻优，其中= ,6。最后我们取得的最优解为(a,b,c)=(11.4980, -1.0430,4.8220),最优值为0.0284，1*10 比上面的情况要差。但是误差可以接受。 我们获得的新的实验对象的变化规律如下: ,1.0430/xy(x),11.4980e，4.8220 (2) 设某热敏电阻的阻值与温度t和内含某稀土元素钍的含量比q的函数关系为非线性最小二乘问题的数学模型为 hgtz(t,x), 1，ht，mq 其中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 试根据背景建立非线性最小二乘问题的数学模型，并用合适的方法求解，得出模型待定参数，分析你的处理过程，你有什么认识, 解: 同上题，我们可以建立问题的最小二乘法模型为 42minf(x),f(x) ,ii,1 hg*10hg*20f(x),,0.126,f(x),,0.219,121，h*10，m*1.01，h*20，m*1.0 hg*10hg*20f(x),,0.076,f(x),,0.126, 321，h*10，m*2.01，h*20，m*2.0 hg*1f(x),,0.18621，h*1，m*0.0 程序: 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); 结果: k，1,6,,f(x),,在实验中，选择的初始点为(0，0，0)，当，停止寻优，其中= 。1*10 7,17(h,g,m),(9.1288，10,1.6121，10,8.9536，10)最终计算得到的最优解为。最优值为 0.0092。 最终的数学模型可以表示为 71.4717，10tz(t,x), 771，9.2188，10t，8.9536，10q 这时我们取第二组数据进行验证，发现在z(20,1),0.1596，和原来的值相差了大概 0.06，偏差有些大了。 造成这种现象的原因可能是因为在建立模型的时候出了问题，换成下面的形式就显得合 理一些。 42minf(x),f(x) ,ii,1 f(x),hg*10,(1，h*10，m*1.0)*0.126,f(x),hg*20,(1，h*20，m*1.0)*0.219,12 f(x),hg*10,(1，h*10，m*2.0)*0.076,f(x),hg*20,(1，h*20，m*2.0)*0.126,32 f(x),hg*1,(1，h*1，m*0.0)*0.1862 这是因为若将约束条件取成上面的除的形式，计算得过程中由于截断误差的原因，会使 得结果变得比较差。但若采取这种和差的形式，可能会好些。我们计算得到的数学模型为(其 ,2,中= ，因为取得再小的话，在自己电脑上运行时，导致迭代的时间太长。不过虽1*10 然精度不是非常高，仍然比上面的情况要好一点。) ,0.0336，tz(t,x),， 1,0.2733，t，0.4509，q 若写成这种形式，效果要好些。由公式计算得到的值如下 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左右，但考虑到这时候的精度 ,2,要比上面小得多= ，所以有理由相信这种做法是好的。(精度取的比较大时，电脑1*10 运行的时间过长。) (3) 许多非数学领域的问题最终也常转化为解线性代数方程组的问题。以系统辨识学科为例，它研究如何通过试验测出系统的输入和输出数据，根据这些数据来建立未知系统的数学模型。设某一被辨识系统的输出y(k)与输入u(k)之间的定量关系为(此模型称为MA 模型)y(k),bu(k)，bu(k,1)，...bu(k,n) 01n 其中诸系数为待定的辨识参数。此模型的参数辨识问题是指:已知输出和输入y(k) 诸信息，确定模型系数。 u(k) 假定在时刻 k=n+1,n+2,L,2n+1 时测量到的数据中不含噪声，那么我们有 y(n，1),bu(n，1)，bu(n)，...bu(1) 01n y(2n，),bu(2n，1)，bu(2n)，...bu(n) 01n 等一系列方程，系统辨识问题转化为求解线性方程组的问题。可以写成下面的形式。如果输入u(k)不是常量，方程有唯一解。 bununuyn(，1)()(1)(，1),,,,,,0,,,,,,bununuyn(，2)(，1)(2)(，2),,,,,,1A,y,A,,,,y, ,,,,,,,, ,,,,,, ,,,,,,bu(2n，1)u(2n)u(n)y(2n，1)n,,,,,, 但是，上述处理基于测量数据中没有噪声的假设，显然与实际情况不符，实际上每次测 量都包含着难以确定的误差, 每个方程需要加上误差，精确解将不存在了。 A,，e,y 为了克服未知噪声的影响，需要增加信息量，比如增加在时刻k=2n+2,2n+3,L,2n+N的测量数据，需要增加信息量，比如增加在时刻k=2n+2,2n+3,L,2n+N 的测量数据，这样我们得到的(1. 6)的系数矩阵A 不是方阵，而是一个“高”的矩阵。出现了一个不相容的矛盾方程组, 这样的线性方程组的解还存在吗,Gauss 提出了最小二乘准则, 用方程的误差 平方和取最小时的为方程的解。 , TTJ,ee,(A,,y)(A,,y) 2,t,t,1/8,t,,tku(t),esin,t，,现在我们假设输入为在取时间间隔的u值，输出为k1 2y(t),10ln(t,1)，,函数在上述时刻的离散值。均为 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 随机正态变量。n=10,N=20, ,,,212 试用最小二乘方法辨识本问题的模型参数，建立本问题的模型。 解: ,然后对下面的无约 问题其实就是一个无约束问题的最优值问题。首先按照题意写出A,y T束问题求最优解。。 ,,(b,b,b...,b)012,n Tmin(A,,y)(A,,y), 其中 u(n，1)u(n)u(1)y(n，1),,,, ,,,,u(n，2)u(n，1)u(2)y(n，2),,,,A,,y,,,,, ,,,, ,,,,u(2n，N)u(2n，N,1)u(n，N)y(2n，N),,,,30，11(30) ,2t2u(k),(esin,t，,)|,y(k),(10ln(t,1)，,)|1t,k,t2t,k,t ,,,其中可由函数randn生成。 12 模型建立完毕。(下面给出代码，代码中描述了整个思想。) 程序: 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