1曲柄连杆结构参数识别.DOC

1、1 曲柄连杆结构参数识别所要识别的参数：R （曲柄半径）1 R=52.5e-3;%曲柄半径2 L=184e-3;%连杆长度3 Ap= 0.0071;%活塞面积4 J=0.3;%曲轴转动惯量5 m=0.4;%单个缸往复质量1.1 简化的转速动力学模型公式推导(把连杆简化成两部分，一部分随活塞上下运动，令一部分随曲轴旋转)-(1)cos1scos11Lrrx根据 lsin其中 为连杆转角， 为曲柄夹角。 。根据三角公式22sin1sin1colr则(1) 变为.2sin1cos(1lrlrx得出 ，带入下式dtrtlcoscs slt(2)rv frrtlt 21 tancoiini 其中 323

2、22 sin14sin1cos;sin1sin gf gfra利用整个曲柄连杆机构瞬时动能相等的原则得：(3)221201vmvJ为换算后的往复和旋转的质量；将(2)带入(2)得到瞬时转动惯量21,m（4）22210sinsinrrmJ然后(3)对 求导得到整个曲柄连杆机构的惯性扭矩T(包括往复惯性扭矩 和旋rT转惯性扭矩 )；eT 2210202001 rmgfrmdJdJd(5)210gfr分为旋转和往复惯性扭矩两部分；往复惯性扭矩 221fTr旋转惯性扭矩 2rme单缸发动机气力扭矩如下 RfApT对于N缸发动机按照发火相位的总和气力扭矩 NkkppfP1惯性扭矩为 N kkkr fgf

3、mT122其中 ， 为转动部分惯量。14kNrotJ根据扭矩平衡关系得到（6）fLrPeTT其中摩擦扭损失压力2105.105.97.0)( Nbartfmep从而得出摩擦扭矩为 fRAptfmeNTf *5).(将上面几个式子带入(6) )7(1212212 fLNkkpNkkkNkk TfPgfrfrmJ 得到二阶转速动力学模型为 )8(1 1211221 xy JfRAxJgxfr fLNkkpNkkk （8）式是一个关于t的非线性二阶微分方程直接求解较为困难，将其转化为关于角度的函数设 ，角域内的 和时域内的 相等，则(9)；tdtt(10)，21dtt令 2x把（10） ， （9）带

4、入（7）得到一阶非齐次线性微分方程（11）QxPd其中,2121NkkfrmJgp212Nkkk fLpfrmJTPRAq1.2 用 matlab 求求微分方程表达式syms o R J L Ap m TL p1 p2 p3 p4 y;f= sin(o)+R*sin(2*o)/(2*(L2-(R*sin(o)2).5);for i=1:4a(i,:)=(subs(f,o,o-(i-1)*pi)2;end;fsum=sum(a);g=diff(f,o);for i=1:4fg(i,:)=subs(f,o,o-(i-1)*pi)*subs(g,o,o-(i-1)*pi);end;fgsum=sum

5、(fg);f1=f;f2=subs(f,o,o-pi);P=2*m*R2*fgsum/(J+m*R2*fsum);Q=2*(Ap*R*(p1+p3)*f1+(p2+p4)*f2)-TL)/(J+m*R2*fsum);dy1=Q- P*y;dy1=(R2*sin(o)2-L2)*(2*TL+y*(-8*m*cos(o)*L2*R2*sin(o)+32*m*cos(o)*R4*sin(o)3-8*m*cos(o)*R4*sin(o)+2*J*cos(o)*R2*sin(o)/(R2*sin(o)2-L2)+(2*R2*cos(o)*sin(o)*(4*m*L2*R2*sin(o)2+J*L2-8*

6、m*R4*sin(o)4+4*m*R4*sin(o)2-J*R2*sin(o)2)/(R2*sin(o)2-L2)2)-2*Ap*R*(p1+p3)*(sin(o)+(R*sin(2*o)/(2*(L2-R2*sin(o)2)(1/2)-(p2+p4)*(sin(o)-(R*sin(2*o)/(2*(L2-R2*sin(o)2)(1/2)/(4*m*L2*R2*sin(o)2+J*L2-8*m*R4*sin(o)4+4*m*R4*sin(o)2-J*R2*sin(o)2)然后令 k(1)=R;k(2)=L;k(3)=Ap;k(4)=J;k(5)=m;k(6)=TLfunction dydo=s

7、im_model(o,y,k,p1,p2,p3,p4)dydo=-(2*k(6)-2*k(3)*k(1)*(p1+p3)*(sin(o)+(k(1)*sin(2*o)/(2*(k(2)2-k(1)2*sin(o)2)(1/2)-(p2+p4)*(sin(o)-(k(1)*sin(2*o)/(2*(k(2)2-k(1)2*sin(o)2)(1/2)/(k(4)-(4*k(1)2*k(5)*sin(o)2*(k(2)2+cos(2*o)*k(1)2)/(k(1)2*sin(o)2-k(2)2)-(2*k(1)2*k(5)*y*(2*k(2)4*sin(2*o)+(5*k(1)4*sin(2*o)/

8、4-k(1)4*sin(4*o)+(k(1)4*sin(6*o)/4-2*k(2)2*k(1)2*sin(2*o)+2*k(2)2*k(1)2*sin(4*o)/(k(1)2*sin(o)2-k(2)2)2*(k(4)-(4*k(1)2*k(5)*sin(o)2*(k(2)2+cos(2*o)*k(1)2)/(k(1)2*sin(o)2-k(2)2);function w=cal_speed_value(k,x)global p1 p2 p3 p4;y0=(76.9280)2;yy=y0;tspan=x;for i= 1:length(tspan)-1t,y=ode45(sim_model,t

9、span(i),tspan(i+1),y0,k,p1(i),p2(i),p3(i),p4(i);y0=y(end,: );yy=yy;y0;endw=yy.0.5;1.3 用求解微分方程的方法识别曲柄连杆结构参数真实值：R=0.0485;%曲柄半径 mL=0.1410;%连杆长度 mAp=0.0058;%活塞面积 m2J=0.23;%转动惯量 kg*m2，带负载或者冷启动时，转动惯量稍大一些，达到 0.38-0.4m=0.33;%单个缸往复质量 kgTL=17.8;%随不同工况变动1.3.1 利用cs#1.xls数据clear;global p1 p2 p3 p4;data=xlsread(c

10、s#1.xls);pdata=data(:,5:8);p11=pdata(203:400,1);p22=pdata(203:400,4);p33=pdata(203:400,3);p44=pdata(203:400,2);p1=(p11+31.4346)*6.895e3;p2=(p22+31.0164)*6.895e3;p3=(p33+29.3030)*6.895e3;p4=(p44+21.1333)*6.895e3;owdata=data(:,2:3);xdata=(owdata(203:400,1)- owdata(203)*pi/180;ydata=owdata(203:400,2);l

11、b=0.02, 0.034, 0.003, 0.15, 0.15, 10;ub=0.1, 0.35, 0.02, 1, 2, 100;k0 = 400, 400, 400, 400, 400,120; options=optimset(TolFun,1e-20,TolX,1e-20,MaxFunEvals,200,Algorithm,trust-region-reflective,Display, iter);kp=lsqcurvefit(cal_speed_value,k0,xdata,ydata,lb,ub,options)kp1=0.0485;0.14;0.0058;0.25;0.15;

12、74%kp= 0.0359 0.0557 0.0087 0.2509 1.0247 83.1993w1= cal_speed_value(kp1,xdata);plot(xdata,ydata,r,xdata,w1,b);1.3.2 利用ne_norm.xls数据clear;global p1 p2 p3 p4;data=xlsread(ne_norm.xls);pdata=data(:,5:8);p11=pdata(363:363+228-1,1);p22=pdata(363:363+228-1,2);p33=pdata(363:363+228-1,3);p44=pdata(363:363+

13、228-1,4);p1=(p11+15.8812+14.5)*6.895e3;p2=(p22+13.8540+14.5)*6.895e3;p3=(p33+12.7006+14.5)*6.895e3;p4=(p44+5.7273+14.5)* 6.895e3;owdata=data(:,2:4);xdata=(owdata(363:363+228-1,1)-owdata(363,1)*pi/180;ydata= owdata(363:363+228-1,2);lb=0.03, 0.034, 0.003, 0.15, 0.15, 10;ub=0.1, 0.35, 0.02, 1, 2, 40;k0

14、 = 0.04, 0.14, 0.06, 0.3, 0.3,30; options=optimset(TolFun,1e-20,TolX,1e-20,MaxFunEvals,100,Algorithm,trust-region-reflective,Display, iter);kp=lsqcurvefit(cal_speed_value,k0,xdata,ydata,lb,ub,options)% kp= 0.0396 0.0683 0.0094 0.6298 0.1503 24.0868w1= cal_speed_value(kp,xdata);plot(xdata,ydata,r,xda

15、ta,w1,b);1.3.3 利用 ne_850.xlsclear;global p1 p2 p3 p4;data=xlsread(ne_850.xls);towdata=data(:,1:4);owdata=data(:,2:4);xdata=(owdata(224:1000,1)-owdata(224,1)*pi/180;ydata= owdata(224:1000,2);pdata=data(:,5:8);p11=pdata(224:1000,1);p22=pdata(224:1000,2);p33=pdata(224:1000,3);p44=pdata(224:1000,4);p1=(

16、p11+16.4886)*6.895e3;p2=(p22+14.6024)*6.895e3;p3=(p33+13.2834)*6.895e3;p4=(p44+4.5723)* 6.895e3;lb=0.03, 0.034, 0.003, 0.15, 0.15, 10;ub=0.1, 0.35, 0.02, 1, 2, 40;k0 = 0.04, 0.14, 0.06, 0.3, 0.3, 30; options=optimset(TolFun,1e-20,TolX,1e-20,MaxFunEvals,100,Algorithm,trust-region-reflective,Display,

17、iter);kp=lsqcurvefit(cal_speed_value,k0,xdata,ydata,lb,ub,options)% kp= 0.0320 0.3468 0.0144 0.5181 0.1500 29.1136w1= cal_speed_value(kp,xdata);plot(xdata,ydata,r,xdata,w1,b);1.4 利用自己改进的扭矩模型进行识别。1.5 final modelTf采用。function y=Tf_model_new211(k,o)y=(k(1)+k(2)*sin(o)+k(3)*sin(2*o)+k(4)*sin(3*o)+k(5)*s

18、in(4*o)+k(6)*sin(o).2+k(7)*sin(2*o).2+k(8)*sin(3*o).2+k(9)*sin(4*o).2;function dydo=sim_model_Tf1(o,y,k,p1,p2,p3,p4)dydo=-(2*(k(6)+k(7)*sin(o)+k(8)*sin(2*o)+k(9)*sin(3*o)+k(10)*sin(4*o)+k(11)*sin(o).2+k(12)*sin(2*o).2+k(13)*sin(3*o).2+k(14)*sin(4*o).2)-2*k(3)*k(1)*(p1+p3)*(sin(o)+(k(1)*sin(2*o)/(2*(

19、k(2)2-k(1)2*sin(o)2)(1/2)-(p2+p4)*(sin(o)-(k(1)*sin(2*o)/(2*(k(2)2-k(1)2*sin(o)2)(1/2)/(k(4)-(4*k(1)2*k(5)*sin(o)2*(k(2)2+cos(2*o)*k(1)2)/(k(1)2*sin(o)2-k(2)2)-(2*k(1)2*k(5)*y*(2*k(2)4*sin(2*o)+(5*k(1)4*sin(2*o)/4-k(1)4*sin(4*o)+(k(1)4*sin(6*o)/4-2*k(2)2*k(1)2*sin(2*o)+2*k(2)2*k(1)2*sin(4*o)/(k(1)2*

20、sin(o)2-k(2)2)2*(k(4)-(4*k(1)2*k(5)*sin(o)2*(k(2)2+cos(2*o)*k(1)2)/(k(1)2*sin(o)2-k(2)2);function w=cal_speed_value_Tf1(k,x)global p1 p2 p3 p4;y0=(71.6349)2;yy=y0;tspan=x;for i= 1:length(tspan)-1t,y=ode45(sim_model_Tf,tspan(i),tspan(i+1),y0,k,p1(i),p2(i),p3(i),p4(i);y0=y(end,: );yy=yy;y0;endw=yy.0.5

21、;clear;global p1 p2 p3 p4;data=xlsread(ne_norm.xls);pdata=data(:,5:8);p11=pdata(363:1000,1);p22=pdata(363:1000,2);p33=pdata(363:1000,3);p44=pdata(363:1000,4);p1=(p11+15.8812)*6.895e3;p2=(p22+13.8540)*6.895e3;p3=(p33+12.7006)*6.895e3;p4=(p44+5.7273)* 6.895e3;owdata=data(:,2:4);xdata=(owdata(363:1000,

22、1)-owdata(363,1)*pi/180;ydata= owdata(363:1000,2);k0=0.0485, 0.14, 0.0058, 0.23, 0.4, -10, 0.1, 0.1, 0.1 -1.4, 1, -0.1, 1 , 1; lb=0.02, 0.034, 0.003, 0, 0, -1e2, -1e3, -1e3, -100, -100, -100, -100, -100, -100;ub=0.1, 0.35, 0.026, 1, 2, 1e2, 1e3, 1e3 100, 100, 100, 100, 1e2, 100;options=optimset(TolF

23、un,1e-20,TolX,1e-20,MaxFunEvals,400,Algorithm,trust-region-reflective,Display, iter);kp,resnorm=lsqcurvefit(cal_speed_value_Tf1,k0,xdata,ydata,lb,ub,options)kp=0.0457 0.1513 0.0064 0.3550 0.1054 4.9133 19.4549 9.1926 2.7378 0.3356 -0.8043 -0.7777 20.9680 1.0000;resnorm= 15.6801;w1= cal_speed_value_T

24、f1(kp,xdata);plot(xdata,ydata,r,xdata,w1,b);2 只识别 J,m 两个参数2.1 把 Tf 当做常数识别（效果不好）function dydo=two_para_model(o,y,k,p1,p2,p3,p4)R=48.5e-3;%曲柄半径 mL=141e-3;%连杆长度 mAp=0.0058;%活塞面积 m2dydo=-(2*k(3)-2*Ap*R*(p1+p3)*(sin(o)+(R*sin(2*o)/(2*(L2-R2*sin(o)2)(1/2)-(p2+p4)*(sin(o)-(R*sin(2*o)/(2*(L2-R2*sin(o)2)(1/2

25、)/(k(1)-(4*R2*k(2)*sin(o)2*(L2+cos(2*o)*R2)/(R2*sin(o)2-L2)-(2*R2*k(2)*y*(2*L4*sin(2*o)+(5*R4*sin(2*o)/4-R4*sin(4*o)+(R4*sin(6*o)/4-2*L2*R2*sin(2*o)+2*L2*R2*sin(4*o)/(R2*sin(o)2-L2)2*(k(1)-(4*R2*k(2)*sin(o)2*(L2+cos(2*o)*R2)/(R2*sin(o)2-L2);function w=cal_speed_value_two(k,x)global p1 p2 p3 p4;y0=(7

26、1.6349)2;yy=y0;tspan=x;for i= 1:length(tspan)-1t,y=ode45(sim_model,tspan(i),tspan(i+1),y0,k,p1(i),p2(i),p3(i),p4(i);y0=y(end,: );yy=yy;y0;endw=yy.0.5;clear;global p1 p2 p3 p4;data=xlsread(ne_norm.xls);pdata=data(:,5:8);p11=pdata(363:363+228-1,1);p22=pdata(363:363+228-1,2);p33=pdata(363:363+228-1,3);

27、p44=pdata(363:363+228-1,4);p1=(p11+15.8812)*6.895e3;p2=(p22+13.8540)*6.895e3;p3=(p33+12.7006)*6.895e3;p4=(p44+5.7273)* 6.895e3;owdata=data(:,2:4);xdata=(owdata(363:363+228-1,1)-owdata(363,1)*pi/180;ydata= owdata(363:363+228-1,2);k0=0.3, 0.3, 30;lb= -1, -2, -100;ub=1, 2, 100;options=optimset(TolFun,1

28、e-20,TolX,1e-20,MaxFunEvals,100,Algorithm,trust-region-reflective,Display, iter);kp=lsqcurvefit(cal_speed_value_two,k0,xdata,ydata,lb,ub,options)w1=cal_speed_value_two(kp,xdata);plot(xdata,ydata,r,xdata,w1,b);方法 2.2 把 Tf 用自己所建立的模型 （效果也不好）Tf采用。function y=Tf_model_new211(k,o)y=(k(3)+k(4)*sin(o)+k(5)*s

29、in(2*o)+k(6)*sin(3*o)+k(7)*sin(4*o)+k(8)*sin(o)2+k(9)*sin(2*o)2+k(10)*sin(3*o)2+k(11)*sin(4*o)2);function dydo=two_para_model_Tf(o,y,k,p1,p2,p3,p4)R=48.5e-3;%曲柄半径 mL=141e-3;%连杆长度 mAp=0.0058;%活塞面积 m2dydo=-(2*(k(3)+k(4)*sin(o)+k(5)*sin(2*o)+k(6)*sin(3*o)+k(7)*sin(4*o)+k(8)*sin(o)2+k(9)*sin(2*o)2+k(10)

30、*sin(3*o)2+k(11)*sin(4*o)2)-2*Ap*R*(p1+p3)*(sin(o)+(R*sin(2*o)/(2*(L2-R2*sin(o)2)(1/2)-(p2+p4)*(sin(o)-(R*sin(2*o)/(2*(L2-R2*sin(o)2)(1/2)/(k(1)-(4*R2*k(2)*sin(o)2*(L2+cos(2*o)*R2)/(R2*sin(o)2-L2)-(2*R2*k(2)*y*(2*L4*sin(2*o)+(5*R4*sin(2*o)/4-R4*sin(4*o)+(R4*sin(6*o)/4-2*L2*R2*sin(2*o)+2*L2*R2*sin(4*

31、o)/(R2*sin(o)2-L2)2*(k(1)-(4*R2*k(2)*sin(o)2*(L2+cos(2*o)*R2)/(R2*sin(o)2-L2);function w=cal_speed_value_two_Tf(k,x)global p1 p2 p3 p4;y0=(71.6349)2;yy=y0;tspan=x;for i= 1:length(tspan)-1t,y=ode45(two_para_model_Tf,tspan(i),tspan(i+1),y0,k,p1(i),p2(i),p3(i),p4(i);y0=y(end,: );yy=yy;y0;endw=yy.0.5;clear;global p1 p2 p3 p4;