储油罐的变位识别.doc

附件一：Matlab程序 程序1：%matlab_1.m%计算值%计算油罐分别在偏移和没偏移时中间部分的每个高度下的实际体积并且绘出图形，同时用计算出来的偏移时的实际体积和附件中实验体积比较close all;clear all;clc;h3=; %代入实测值，高度v3= %代入实测值，体积v5=;v5=v5+262;h5=;a=0.89;b=0.6;l=2.45;a1=4.1/180*pi;l1=0.825;h1=147:1:1171; %中间部分的H的范围h1=h1/1000 %小数号h=h1-tan(a1)*l1; % h代入forfor i=1:1:1025 %有偏移，变位后，计算值v(i)=a*b*l*( pi/2+(h(i)/b-1)*sqrt(2*h(i)/b-(h(i)/b)2)+asin(h(i)/b-1)-(1/12*(h(i)/b-1)/sqrt(1-(h(i)/b-1)2)*(tan(a1)2*l3*a/b;end;v=v*1000; % v :1000 for i=1:1:1025 % a=0没变位，计算值v1(i)=a*b*l*(pi/2+(h1(i)/b-1)*sqrt(2*h1(i)/b-(h1(i)/b)2)+asin(h1(i)/b-1) ;end;v1=v1*1000; %1000h1=h1*1000; % 1000 h3=h3/1000; h4=h3-tan(a1)*l1;for i=1:1:53 %把附件中实验高度代入偏移的函数中求出体积c_v(i)=a*b*l*(pi/2+(h4(i)/b-1)*sqrt(2*h4(i)/b-(h4(i)/b)2)+asin(h4(i)/b-1)-(1/12*(h4(i)/b-1)/sqrt(1-(h4(i)/b-1)2)*(tan(a1)2*l3*a/b;end;c_v=c_v*1000 ; %理论计算值，53点c_v=c_v;v3; %实测值dt_v=c_v-v3; %计算值与实测值的误差for i=1:1:53 wc(i)=dt_v(i)/c_v(i); %计算值与实测值的误差百分比 endwc=wc*100;ave1=sum(wc)/53;h5=h5/1000;for i=1:1:78 % a=0没变位，计算值v6(i)=a*b*l*(pi/2+(h5(i)/b-1)*sqrt(2*h5(i)/b-(h5(i)/b)2)+asin(h5(i)/b-1) ;end;v6=v6*1000;v6h5=h5*1000; dt_v5=v6-v5; %计算值与实测值的误差dt_v5;for i=1:1:78 wc5(i)=dt_v5(i)/v6(i); %计算值与实测值的误差百分比endwc5=wc5*100ave2=sum(wc5)/78;ave2v_mid_low=v(1)v_mid_up=v(1025)figureplot(h5,wc5,b*); grid on;h3=h3*1000;figureplot(h3,wc,bo); grid on;figureplot(h1,v,r); %有位差，理论值 grid on;hold on;plot(h3,v3,k); %画实测值，53点grid on;hold on;plot(h1,v1,g); %无位差值 hold onplot(h5,v5,b); %无位差值 hold off;程序2：%matlab_2.m%分析变位对罐容表的影响clear all;clc;a=0.89;b=0.6;l=2.45;a1=4.1/180*pi;l1=0.825;h1=147:1:1171;h1=h1/1000;h=h1-tan(a1)*l1;for i=1:1:1025 %变位后，理论计算值v(i)=a*b*l*(pi/2+(h(i)/b-1)*sqrt(2*h(i)/b-(h(i)/b)2)+asin(h(i)/b-1)-(1/12*(h(i)/b-1)/sqrt(1-(h(i)/b-1)2)*(tan(a1)2*l3;end;v=v*1000;v;for i=1:1:1025 %不变位，理论计算值v1(i)=a*b*l*(pi/2+(h1(i)/b-1)*sqrt(2*h1(i)/b-(h1(i)/b)2)+asin(h1(i)/b-1) ;end;v1=v1*1000;h1=h1*1000;v1;v2=v1-v; %同一高度下有变位和没变位的体积差v2ave=sum(v2)/907;avefigureplot(h1,v,r); %有位差，理论值grid on;hold on;plot(h1,v1,g); %无位差值 hold off;figureplot(v2,h1,g)grid on;程序3：%matlab_3.m%模型建立后a=0和a=4.1的曲线close all;clear all;clc;a=0.89;b=0.6;l=2.45;a1=4.1/180*pi;l1=0.825;h9=1:1:1200;h9=h9/1000;for i=1:1:1200 %a=0没变位，计算值v9(i)=a*b*l*(pi/2+(h9(i)/b-1)*sqrt(2*h9(i)/b-(h9(i)/b)2)+asin(h9(i)/b-1) ;end;h9=h9*1000;v9=v9*1000;hx=2.05*tan(a1)hy=1.2-0.4*tan(a1)hx=hx*1000;hy=hy*1000;hz=hy-hx+1h1=hx:1:hy; %中间部分的H的范围h1=h1/1000 ; %小数号h=h1-tan(a1)*l1; % h代入 forfor i=1:1:hz %有偏移，变位后，计算值v(i)=a*b*l*(pi/2+(h(i)/b-1)*sqrt(2*h(i)/b-(h(i)/b)2)+asin(h(i)/b-1)-(1/12*(h(i)/b-1)/sqrt(1-(h(i)/b-1)2)*(tan(a1)2*l3*a/b;end;v=v*1000; % v :1000 h1=h1*1000; % 1000 h3=0:1:hx; %油罐底面部分H的范围是0147mmh3=h3/1000; l2=0.4+h3/tan(a1);l3=l2/2-0.4;h2=h3-tan(a1)*l3;for i=1:1:(hx+1) %每个高度下的实际体积 v3(i)=a*b*l2(i)*(pi/2+(h2(i)/b-1)*sqrt(2*h2(i)/b-(h2(i)/b)2)+asin(h2(i)/b-1)-(1/12*(h2(i)/b-1)/sqrt(1-(h2(i)/b-1)2)*(tan(a1)2*l2(i)3*a/b;end;v3=v3*1000; h3=h3*1000; hm=1.65*tan(a1)*1000h4=hx:-1:(hm+1);h4=h4/1000; %小数号l4=0.4+h4/tan(a1);l5=l4/2-0.4;h6=h4-tan(a1)*l5;hn=hx-hmfor i=1:1:hn %变位后，理论计算值，1200点v7(i)=a*b*l4(i)*(pi/2+(h6(i)/b-1)*sqrt(2*h6(i)/b-(h6(i)/b)2)+asin(h6(i)/b-1)-(1/12*(h6(i)/b-1)/sqrt(1-(h6(i)/b-1)2)*(tan(a1)2*l4(i)3)*a/b;end;v8=pi*a*b*2.45-v7; %算油罐偏移时上面部分的每个高度下的实际体积v8=v8*1000; % v :1000 h4=h4*1000; % 1000 h4=1200-h4+118; %上面部分的高度H的范围是11711200mmfigureplot(h9,v9,b); %有位差，理论值hold on;plot(h1,v,r); %有位差，理论值hold on;plot(h3,v3,b); %绘出曲线hold on;plot(h4,v8,b); grid on; 程序4：%matlab_4.m%非线性回归，函数，function v3 function v3=my_a_b(beta,x) h11=x(:,1);h12=x(:,2);r=1.5;rr=1.625;m=0.625;l=8;l0=2;n=3000;h11=r-(r-h11)*cos(beta(2)-l0*tan(beta(1); h1=h11+tan(beta(1)*l/2; h2=h11-tan(beta(1)*l/2;i1=h1*1000;i2=h2*1000;i3=h11;vhat1=r2*l*(pi/2+(i3/r-1).*sqrt(2*i3/r-(i3/r).2)+asin(i3/r-1) -. %体积v函数 1/12*(i3/r-1)./sqrt(1-(i3/r-1).2)*(tan(beta(1)2*l3 +. r/(3*n3)*(2*i1-n).*(3*rr2*n2-r2*(2*i1-n).2).*asin(2*r*sqrt(i1.*(n-i1)./sqrt(n2*rr2-r2*(2*i1-n).2)+2*rr3/3*(pi/2+atan(m*(2*i1-n)./(2*rr*sqrt(i1.*(n-i1)-m/3*(2*rr2+r2).*(pi/2+asin(2*i1-n)./n)-4*m*r2/(3*n2)*(2*i1-n).*sqrt(i1.*(n-i1)+. r/(3*n3)*(2*i2-n).*(3*rr2*n2-r2*(2*i2-n).2).*asin(2*r*sqrt(i2.*(n-i2)./sqrt(n2*rr2-r2*(2*i2-n).2)+2*rr3/3*(pi/2+atan(m*(2*i2-n)./(2*rr*sqrt(i2.*(n-i2)-m/3*(2*rr2+r2).*(pi/2+asin(2*i2-n)./n)-4*m*r2/(3*n2)*(2*i2-n).*sqrt(i2.*(n-i2);h12=r-(r-h12)*cos(beta(2)-l0*tan(beta(1); h1=h12+tan(beta(1)*l/2; h2=h12-tan(beta(1)*l/2;i1=h1*1000;i2=h2*1000;i3=h12;vhat2=r2*l*(pi/2+(i3/r-1).*sqrt(2*i3/r-(i3/r).2)+asin(i3/r-1) -. %体积v函数 1/12*(i3/r-1)./sqrt(1-(i3/r-1).2)*(tan(beta(1)2*l3 +. r/(3*n3)*(2*i1-n).*(3*rr2*n2-r2*(2*i1-n).2).*asin(2*r*sqrt(i1.*(n-i1)./sqrt(n2*rr2-r2*(2*i1-n).2)+2*rr3/3*(pi/2+atan(m*(2*i1-n)./(2*rr*sqrt(i1.*(n-i1)-m/3*(2*rr2+r2).*(pi/2+asin(2*i1-n)./n)-4*m*r2/(3*n2)*(2*i1-n).*sqrt(i1.*(n-i1)+. r/(3*n3)*(2*i2-n).*(3*rr2*n2-r2*(2*i2-n).2).*asin(2*r*sqrt(i2.*(n-i2)./sqrt(n2*rr2-r2*(2*i2-n).2)+2*rr3/3*(pi/2+atan(m*(2*i2-n)./(2*rr*sqrt(i2.*(n-i2)-m/3*(2*rr2+r2).*(pi/2+asin(2*i2-n)./n)-4*m*r2/(3*n2)*(2*i2-n).*sqrt(i2.*(n-i2);v3=vhat1-vhat2;v3=real(v3);程序5：%matlab_5.m%非线性回归，主函数，求alpha,betaclear all;close all;clc;%format longh11=xlsread(D:data_second.xls, 1, E2:E302); %横坐标h11read_h=1h12=xlsread(D:data_second.xls, 1, E3:E303); %横坐标h12v3=xlsread(D:data_second.xls, 1, d3:d303); %纵坐标v3h11=h11/1000;h12=h12/1000;v3=v3/1000; x(:,1)=h11; % mx(:,2)=h12;y=v3; % m3beta=3,5.04; %回归系数的初值,任意取的 beta=beta/180*pi; %pi转换betahat,r,j=nlinfit(x,y,my_a_b,beta); %r,j是下面命令用的信息 betaci=nlparci(betahat,r,jacobian,j); r; % 残差%plot(x,r,+) % 残差图betaa=betahat,betaci %回归系数及其置信区间yhat,delta=nlpredci(my_a_b,x,betahat,r,jacobian,j)%plot(x,y,+,x,yhat,r) % 预测及作图%y的预测值及其置信区间的半径，置信区间为yhatdeltanlintool(x,y,my_a_b,beta); %用nlintool得到一个交互式画面，左下方%的Export可向工作区传送数据，如剩余标准差等%betaa 回归系数及其置信区间%yhatdelta y的预测值及其置信区间的半径，置信区间为 程序6：%matlab_6.m%求变位后，罐容表的标定值clear all; clc;h=0:0.1:3; %求取03m,步进为10cm的高对应体积a1=2.1257;b1=4.6868;beta(1)=a1; %角度1beta(2)=b1; %角度2beta=beta/180*pi;r=1.5;rr=1.625;m=0.625;l=8;l0=2;n=3000;h=r-(r-h)*cos(beta(2)-l0*tan(beta(1); %beta(2)h1=h+tan(beta(1)*l/2; h2=h-tan(beta(1)*l/2;i1=h1*1000;i2=h2*1000;i3=h;vhat=r2*l*(pi/2+(i3/r-1).*sqrt(2*i3/r-(i3/r).2)+asin(i3/r-1) -. %体积v函数 1/12*(i3/r-1)./sqrt(1-(i3/r-1).2)*(tan(beta(1)2*l3 +. r/(3*n3)*(2*i1-n).*(3*rr2*n2-r2*(2*i1-n).2).*asin(2*r*sqrt(i1.*(n-i1)./sqrt(n2*rr2-r2*(2*i1-n).2)+2*rr3/3*(pi/2+atan(m*(2*i1-n)./(2*rr*sqrt(i1.*(n-i1)-m/3*(2*rr2+r2).*(pi/2+asin(2*i1-n)./n)-4*m*r2/(3*n2)*(2*i1-n).*sqrt(i1.*(n-i1)+. r/(3*n3)*(2*i2-n).*(3*rr2*n2-r2*(2*i2-n).2).*asin(2*r*sqrt(i2.*(n-i2)./sqrt(n2*rr2-r2*(2*i2-n).2)+2*rr3/3*(pi/2+atan(m*(2*i2-n)./(2*rr*sqrt(i2.*(n-i2)-m/3*(2*rr2+r2).*(pi/2+asin(2*i2-n)./n)-4*m*r2/(3*n2)*(2*i2-n).*sqrt(i2.*(n-i2);vhat=real(vhat);vhat=vhat程序7:%matlab_7.m%目标函数曲线与实测数据曲线比较%my_a_b.m： clear all;clc;h=; %代入数据v=; %代入数据h=h/1000;v=v/1000;a1=0; %令角度为0b1=0; %令角度为0beta(1)=a1;beta(2)=b1;beta=beta/180*pi;r=1.5;rr=1.625;m=0.625;l=8;l0=2;n=3000;h=r-(r-h)*cos(beta(2)-l0*tan(beta(1); % beta(2)h1=h+tan(beta(1)