工程数学作业第十四次方健(12.29).docx

姓名：方健 学号：652081701073 专业：化学工程 班级：化工1503问题：在一个等温间歇反应器中进行复杂反应，描述该反应系统的七个微分方程共有5个参数k1，k2，k3，k4和k5，初试状态x（0）=0.1883 0.2507 0.0467 0.0899 0.1804 0.1394 0.1046T,试根据下表的数据，利用MATLAB最优化方法函数估计出k参数，表中的数据存于文件KineticsData1.m中% 动力学数据:% t y1 y2 y3 y4 ExpData = . 0 0.18830.0899 0.18040.1394 0.0100 0.2047 0.0866 0.1729 0.1297 0.0200 0.2181 0.0856 0.1680 0.1205 0.0300 0.2291 0.0863 0.1647 0.1123 0.0400 0.2382 0.0878 0.1623 0.1053 0.0500 0.2459 0.0899 0.1604 0.0995 0.0600 0.2523 0.0921 0.1588 0.0948 0.0700 0.2576 0.0945 0.1574 0.0911 0.0800 0.2622 0.0968 0.1561 0.0882 0.0900 0.2660 0.0989 0.1548 0.0859 0.1000 0.2692 0.1010 0.1537 0.0842 0.1100 0.2719 0.1028 0.1525 0.0830 0.1200 0.2742 0.1045 0.1515 0.0821 0.1300 0.2761 0.1060 0.1505 0.0814 0.1400 0.2777 0.1074 0.1495 0.0810 0.1500 0.2790 0.1086 0.1487 0.0807 0.1600 0.2801 0.1096 0.1479 0.0805 0.1700 0.2811 0.1106 0.1471 0.0803 0.1800 0.2819 0.1114 0.1465 0.0803 0.1900 0.2825 0.1121 0.1458 0.0803 0.2000 0.2830 0.1127 0.1453 0.0803x = lsqnonlin(fun,x0)x = lsqnonlin(fun,x0,lb,ub)x = lsqnonlin(fun,x0,lb,ub,options)x = lsqnonlin(problem) x,resnorm,residual,exitflag,output,lambda,jacobian = lsqnonlin(.)x = lsqcurvefit(fun,x0,xdata,ydata)x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub)x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)x = lsqcurvefit(problem)x,resnorm,residual,exitflag,output,lambda,jacobian = lsqcurvefit(.)根据问题首先进行方程参数的估计编写KineticsEst5,m文件如下：function KineticsEst5% 动力学ODE方程模型的参数估计% The variables y here are y(1)=x1, y(2)=x4, y(3)=x5,y(4)=x6 . clear allclck0 = 0.5 0.5 0.5 0.5 0.5; % 参数初值lb = 0 0 0 0 0; % 参数下限ub = +inf +inf +inf +inf +inf; % 参数上限x0 = 0.1883 0.2507 0.0467 0.0899 0.1804 0.1394 0.1046;KineticsData1;yexp = ExpData(:,2:5); % yexp: 实验数据x1x4x5x6% 使用函数fmincon()进行参数估计k,fval,flag = fmincon(ObjFunc4Fmincon,k0,lb,ub,x0,yexp);fprintf(n使用函数fmincon()估计得到的参数值为:n)fprintf(tk1 = %.4fn,k(1)fprintf(tk2 = %.4fn,k(2)fprintf(tk3 = %.4fn,k(3)fprintf(tk4 = %.4fn,k(4)fprintf(tk5 = %.4fn,k(5)fprintf( The sum of the squares is: %.1enn,fval)k_fmincon = k;% 使用函数lsqnonlin()进行参数估计k,resnorm,residual,exitflag,output,lambda,jacobian = . lsqnonlin(ObjFunc4LNL,k0,lb,ub,x0,yexp); kci = nlparci(k,residual,jacobian)fprintf(nn使用函数lsqnonlin()估计得到的参数值为:n)output% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计k0 = k_fmincon;k,resnorm,residual,exitflag,output,lambda,jacobian = . lsqnonlin(ObjFunc4LNL,k0,lb,ub,x0,yexp); kci = nlparci(k,residual,jacobian)fprintf(nn以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:n)output% -function f = ObjFunc4Fmincon(k,x0,yexp)tspan = 0.00 : 0.01 : 0.20;t x = ode45(KineticEqs,tspan,x0,k); y(:,1) = x(:,1);y(:,2:4) = x(:,4:6);f = sum(y(:,1)-yexp(:,1).2) + sum(y(:,2)-yexp(:,2).2) . + sum(y(:,3)-yexp(:,3).2) + sum(y(:,4)-yexp(:,4).2);% -function f = ObjFunc4LNL(k,x0,yexp)tspan = 0.00 : 0.01 : 0.20; t x = ode45(KineticEqs,tspan,x0,k); y(:,1) = x(:,1);y(:,2:4) = x(:,4:6);f1 = y(:,1) - yexp(:,1); f2 = y(:,2) - yexp(:,2);f3 = y(:,3) - yexp(:,3);f4 = y(:,4) - yexp(:,4);f = f1; f2; f3; f4;% -function dxdt = KineticEqs(t,x,k)q = 8.75 + k(5);dxdt = . ( k(5)-q*x(1)- k(1)*x(1)*x(2)-k(4)*x(1)*x(6)*sqrt(0.9) ) ( 7.0-q*x(2) - k(1)*x(1)*x(2)-2*k(2)*x(2)*x(3) ) ( 1.75 -q*x(3) - k(2)*x(2)*x(3) ) ( -q*x(4) + 2*k(1)*x(1)*x(2)-k(3)*x(4)*x(5) ) ( -q*x(5) + 3*k(2)*x(2)*x(3)-k(3)*x(4)*x(5) ) ( -q*x(6) + 2*k(3)*x(4)*x(5)-k(4)*x(1)*x(6)*sqrt(0.9) ) ( -q*x(7) + 2*k(4)*x(1)*x(6)*sqrt(0.9) );在MATLAB中运行结果如下：ExpData = 0 0.1883 0.0899 0.1804 0.1394 0.0100 0.2047 0.0866 0.1729 0.1297 0.0200 0.2181 0.0856 0.1680 0.1205 0.0300 0.2291 0.0863 0.1647 0.1123 0.0400 0.2382 0.0878 0.1623 0.1053 0.0500 0.2459 0.0899 0.1604 0.0995 0.0600 0.2523 0.0921 0.1588 0.0948 0.0700 0.2576 0.0945 0.1574 0.0911 0.0800 0.2622 0.0968 0.1561 0.0882 0.0900 0.2660 0.0989 0.1548 0.0859 0.1000 0.2692 0.1010 0.1537 0.0842 0.1100 0.2719 0.1028 0.1525 0.0830 0.1200 0.2742 0.1045 0.1515 0.0821 0.1300 0.2761 0.1060 0.1505 0.0814 0.1400 0.2777 0.1074 0.1495 0.0810 0.1500 0.2790 0.1086 0.1487 0.0807 0.1600 0.2801 0.1096 0.1479 0.0805 0.1700 0.2811 0.1106 0.1471 0.0803 0.1800 0.2819 0.1114 0.1465 0.0803 0.1900 0.2825 0.1121 0.1458 0.08030.2000 0.2830 0.1127 0.1453 0.0803使用函数fmincon()估计得到的参数值为:k1 = 17.5835k2 = 72.7641k3 = 51.1175k4 = 22.9402k5 = 5.9949 The sum of the squares is: 6.4e-07k = 17.6085 73.0645 51.3262 23.0250 6.0012ci = 17.5954 17.6215 72.9804 73.1485 51.2783 51.3742 22.9682 23.0817 5.9989 6.0036使用函数lsqnonlin()估计得到的参数值为:output = firstorderopt: 2.6911e-07 iterations: 9 funcCount: 60 cgiterations: 0 algorithm: trust-region-reflective message: 1x425 chark = 17.6063 73.0504 51.3184 23.0156 6.0009ci = 17.5933 17.6194 72.9662 73.1345 51.2704 51.3663 22.9588 23.0724 5.9985 6.0032以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:output = firstorderopt: 8.3868e-07 iterations: 1 funcCount: 12 cgiterations: 0 algorithm: trust-region-reflective message: 1x425 char在微分法的情况下：function KineticsEst1_int% 动力学参数辨识: 用积分法进行反应速率分析得到速率常数k和反应级数n% Analysis of kinetic rate data by using the integral method% Reaction of the type - rate = kCAorder% order - reaction order% rate - reaction rate vector% CA - concentration vector for reactant A% T - vector of reaction time% N - number of data points% k- reacion rate constantclear allclcglobal CAmt = 0 20 40 60 120 180 300;CAm = 10 8 6 5 3 2 1;% 非线性拟合beta0 = 0.0053 1.39;tspan = 0 20 40 60 120 180 300;CA0 = 10;beta,resnorm,resid,exitflag,output,lambda,jacobian = . lsqnonlin(OptObjFunc,beta0,tspan,CA0,CAm)ci = nlparci(beta,resid,jacobian)% 拟合效果图(实验与拟合的比较)t4plot CA4plot = ode45(KineticsEqs,tspan(1) tspan(end),CA0,beta);plot(t,CAm,bo,t4plot,CA4plot,k-)legend(Exp,Model)xlabel(时间t, s)ylabel(浓度C_A, mol/L)% 残差关于拟合值的残差图t CAc = ode45(KineticsEqs,tspan,CA0,beta);figureplot(CAc,resid,*)xlabel(浓度拟合值（mol/L）)ylabel(残差R (mol/L)refline(0,0)% 参数辨识结果fprintf(Estimated Parameters:n)fprintf(tk = %.4f %.4fn,beta(1),ci(1,2)-beta(1)fprintf(tn = %.2f %.2fn,beta(2),ci(2,2)-beta(2)% -function f = OptObjFunc(beta,tspan,CA0,CAm)t CAc = ode45(KineticsEqs,tspan,CA0,beta);f = CAc - CAm;% -function dCAdt = Kine