学士学位论文--连串平行及连串反应的等温优化.doc

毕 业 论 文

题 目：连串—平行及连串反应的等温优化
学 院： 物理与化学学院
年级、专业： 2009级、化学
学 生： 黄芮
学 号： 312009070301236
指导教师： 彭昌荣
完成日期： 2013年 5月 20日

2

西华大学本科毕业论文
目 录

摘 要 3
Abstract 4
前 言 5
1 等温优化计算实例 5
1.1 等温优化结果 7
2 非等温优化计算实例 12
2.1 对该反应作等温搜索的计算结果 15
3 结果与讨论 16
总结与体会 18
谢 辞 19
参考文献 20
附MATLAB程序 21

摘 要
在精细化工生产中，以复杂反应系统为特征的、能生成不希望的最终产物的连串-平行反应或连串反应是俯拾皆是的，而这一类反应的目标又常常是希望获得高质量和高纯度的产品，根据反应的特点和约束条件，进行操作条件的优化是必要的。本文首先对一个连串-平行反应实例的操作条件进行了等温优化并编写了相应的Matlab程序；其次对一个经典的连串反应的温度控制问题，采用各间隔反应时间内等温的方法进行了目标产物浓度最大问题的优化并编写了相应的Matlab程序；第二个计算实例中提供的优化算法比庞特里亚金极大值原理要简单的多，尽管其优化的温度分布是次优的，但是，可以通过改变程序中间隔反应时间的数目来逼近最优。
关键词：连串-平行反应 连串反应 等温反应时间间隔 庞特里亚金极大值原理

Abstract
Consecutive-parallel or consecutive reaction is still widely used in fine chemical production which is characterized by quite complex reaction system which can produce undesirable end products. As the aim of the fine chemical industry is to produce high quality and purity products,It is essential to optimize operating conditions，taking into account constrains and reation’s characteristics. In this work, the first part of this paper dealed with the optimal isothermal operating conditions of a consecutive-parallel and the corresponding Matlab program for this type of reaction was programmed. Then the optimal temperature profile of a classical consecutive, taking into account the mximum concentration of objective product, and the corresponding Matlab program for this type of reaction were carried out by the several isothermal subintervals of reaction time on the basis of the characteristics of the studied reaction. The second computing example’s method proposed by myself compared with Pontryagin maximum principle was relatively simple. By Increasing the number of isothermal subintervals of reaction time, I’m sure the optimal results can approach to Pontryagin maximum principle’s computing results.
Keywords: Consecutive-parallel reaction; Consecutive reaction; Isothermal subintervals of reaction time;Pontryagin maximum principle.

前 言
连串(平行)反应是反应产物可以进一步反应生成其他产物的反应, 它是化学工业中常见的一类复杂反应[1,2]。如氯化苯合成、烷基苯合成、以硝基苯为原料的对氨基苯酚的电解合成等均属于此类反应[3]。在实际工艺生产中, 未反应的原料总是要回收利用。 连串反应工艺优化的基本目标是使消耗的原料尽可能多的得到目的产物根据优化的具体目标不同, 有不同的优化处理方法，连串反应优化问题有[4,5]：等温优化，即给定反应时间，优化反应温度。非等温优化[6]、最佳进料比优化以及成本最低优化法[7]。本文讨论给定反应时间，要求目的产物浓度最大时优化反应温度，即求出反应温度与反应时间的对应（数值）关系。
1 等温优化计算实例
采用氨() 与环氧丙烷(PO) 为原料, 以高氨环比(NH3/PO) 先选择性制备一异丙醇胺(MIPA), 再由MIPA 与PO 合成二异丙醇胺( DIPA) 是最近发展起来的新工艺[8] , 它与传统的生产工艺相比能有效地提高DIPA 的生产选择性, 降低副产三异丙醇胺( TIPA) 的生成。蒋旭峰, 曾崇余, 任晓乾等研究了由一异丙醇胺合成二异丙醇胺的反应规律[9,10] , 在文献[9,10]的基础上对该反应进行了动力学方程方面初步研究, 为工业化放大提供基础数据[11]，得出MIPA与PO的反应为一连串-平行反应（consecutive-parallel or consecutive-competitive reaction）， 如下所示：
MIPA+ PODIPA （1）
DIPA+ POTIPA （2）
体系中各组分反应级数均为一级，体系中各组份的反应速度可表示为：
简记为： （3）
简记为： （4）
简记为： （5）
其中：
初始条件： （6）
（7）
单位为：
（8）
单位为：
适用范围：MIPA 与PO 的摩尔比在1：1~ 1：3, 反应温度在30~70℃ 的实验数据作为拟合样本的，经检验，实验值与模型计算值能较好地吻合， 其相对偏差小于6%， 说明该动力学方程是可靠的。
其中，PO：表示环氧丙烷；MIPA：表示一异丙醇胺；DIPA：表示二异丙醇胺；TIPA：表示三异丙醇胺。
ｃ：表示体积摩尔浓度， 。
A1：表示反应（1） 的指前因子，。
A2：表示反应（2） 的指前因子，。
：表示反应（1） 的活化能， 。
：表示反应（2） 的活化能，。
R：表示通用气体常数， 8. 314。
T：表示温度，K。
这样得状态方程（浓度或质量平衡方程）：
（9）
1.1 等温优化结果
优化参数为给定反应物初始浓度条件下的温度（温度区间[20 70]℃）和反应时间，在反应物初始浓度条件为:
的条件下，使得目的产物B的浓度（以下以表示）最大的最佳反应温度。
式（9）是一个一阶常微分方程的初值问题，在计算机普及的今天，有好几种数学软件可以求解此类问题[12]，本文应用MATLAB来解此类问题，编程和作图都相对简单。
若反应温度为30℃，反应时间500min，通过计算得的分布见图1，反应时间250min，才接近最大值，以后反而下降，见图2。

图1 反应温度30℃，反应时间500min时系统中各组分浓度-时间分布

图2 反应温度30℃，反应时间500min时浓度-时间分布
若反应温度为50℃，反应时间40min，的分布见图3，反应时间30min，才接近最大值，以后反而下降，见图4。可见反应温度的提高使得反应时间大大缩短，而都接近1.4

图3 反应温度50℃，反应时间40min时系统中各组分浓度-时间分布

图4 反应温度50℃，反应时间40min时浓度-时间分布
若反应温度为60℃，反应时间40min，的分布见图5，反应时间10min，接近最大值，以后反而下降。可见反应温度的提高使得反应时间大大缩短，而都接近1.4，所以，优化温度可以取60℃，达到的反应时间不会超过10min。

图5 反应温度60℃，反应时间40min时浓度-时间分布
若取反应时间为10min，可以计算各反应温度下的最大值，与反应温度T(℃)的关系即的分布见图6，由图6可见反应时间10min，温度大于60℃以后提高温度对增加效果甚微。于是该初始浓度条件下的优化温度可以取60℃，反应时间取10min，而都接近1.4，所以，优化温度可以取60℃，达到的反应时间不会超过10min。反应温度60℃，反应时间10min时系统各组分浓度分布见图7，最大值。其余初始浓度条件下的优化温度和反应时间的优化仿此。

图6 反应时间10min时各反应温度下最大值 -时间分布

图7 反应温度60℃，反应时间10min时系统各组分浓度分布

图8 反应温度60℃，反应时间10min时浓度-时间分布
2 非等温优化计算实例
下面以文献[13,14,15]研究了的反应，，其中，，温度范围：,初始条件：，，。目标函数是研究在给定的间歇反应时间时，求出一个温度分布使得B的浓度最大，即。不同的是本文采用自己提出的方法：时间区间为tspan=[0.00:1/n:1.00]，时间间隔 ，假设初始温度分布为Temp0 =398:-(398-298)/n:298，并且控制过程中假定温度是由高到低的，因为第一反应的活化能低于后一个反应的活化能，所以前期反应温度高后期温度低会相对有利，在每一个时间内假设是等温反应，赋予一个反应温度，在[0 dt]时间范围积分得到各个组分浓度分布，下一个时间内赋予另一个反应温度，并且初始组成为上一个时间结束时的组成，在[0 dt]时间范围积分得到各个组分浓度分布，一次循环，直到tspan的终点，然后寻找的温度分布即为最优温度分布。计算20个时间间隔得到的使得最大的温度分布及组成分布数据见表。若改变时间间隔数目，温度分布会不一样，这取决于工艺对控制的要求；若温度被由低到高来控制（改变温度约束矩阵A中的-1为1），会得到与下面的等温搜索差不多相同的结果；若温度控制可以忽高忽低，又会有不同的温度分布；一句话温度控制策略不同会有不同的温度分布。
表1 20个时间间隔计算得到的使得最大的温度分布及组成分布数据
Time/s
T/K

0.00
363.3
1
0
0
0.00-0.05
363.3
0.8295
0.1676
0.0029
0.05-0.10
351.4
0.7308
0.2618
0.0074
0.10-0.15
345.8
0.6608
0.3270
0.0122
0.15-0.20
343.1
0.6059
0.3767
0.0173
0.20-0.25
341.6
0.5608
0.4165
0.0228
0.25-0.30
340.8
0.5225
0.4490
0.0285
0.30-0.35
340.1
0.4897
0.4759
0.0344
0.35-0.40
339.2
0.4612
0.4984
0.0404
0.40-0.45
338.1
0.4364
0.5172
0.0464
0.45-0.50
336.6
0.4149
0.5330
0.0521
0.50-0.55
334.7
0.3961
0.5463
0.0576
0.55-0.60
332.5
0.3797
0.5576
0.0627
0.60-0.65
329.9
0.3655
0.5672
0.0672
0.65-0.70
326.9
0.3532
0.5755
0.0713
0.70-0.75
323.6
0.3425
0.5827
0.0748
0.75-0.80
320.0
0.3332
0.5890
0.0778
0.80-0.85
316.2
0.3253
0.5945
0.0803
0.85-0.90
312.1
0.3184
0.5993
0.0823
0.90-0.95
307.8
0.3125
0.6036
0.0840
0.95-1.00
303.3
0.3074
0.6073
0.0853
1.00
298.7
0.3031
0.6107
0.0863

图9 温度分布图（阶梯图,20个时间间隔）

图10 温度分布图

图11 ()随反应时间的分布
2.1 对该反应作等温搜索的计算结果
等温反应搜索到的见表和图，可见若为等温反应的话，达到最大的温度区间在330~340℃，335℃是比较合适的，样条插值法内插得到最佳温度为335.3℃，与最优温度分布计算的到的=0.6107差别很小，而且对于温度的控制来说难度也相对较小，而非等温反应对控制系统要求是比较高的，但是能够将二者计算出来做一对比也是很有意义的，尽管非等温反应控制比较难，也可能对于其他情况用非等温是必要的。一般是最优温度分布的结果略大于等温优化的结果，在二者差异不大时，宁可采用相对简单的控制策略是明智的选择。
表2 等温搜索数据
T/K
398
393
388
383
378
373
368

0.4579
0.4699
0.4823
0.4946
0.5067
0.5192
0.5320
T/K
363
358
353
348
343
338
333

0.5450
0.5583
0.5717
0.5854
0.5982
0.6050
0.6053
T/K
328
323
318
313
308
303
298

0.5994
0.5879
0.5714
0.5504
0.5257
0.4977
0.4671

图12 等温反应时的分布
3 结果与讨论
（1）对连串-平行反应实例作了等温优化，并得到了需要的结果。
（2）对连串反应实例作了非等温优化，并采用了自己提出的次优办法，得到了需要的结果。该方法比应用Pontryagin最大值原理更简单有效，况且可以改变等温反应时间间隔的数目来进一步逼近Pontryagin最大值原理的结果，计算结果表明，时间间隔的数目的增加到一定程度后计算所得的温度分布对于增加目标产物的浓度意义并不大，而于实际控制来说，难度会急剧增加，这可能是不经济的。同时对该反应作了等温优化并与非等温优化结果作了比较，二者相差很小（非等温优化结果，等温优化结果），对于非等温反应来说温度的控制难度较大，对等温反应温度的控制难度相对较小。
（3）对于同一反应能够计算出等温与非等温优化的结果并进行比较，以便做出温度控制的选择是很有意义的，那样我们对反应会做到“心中有数”。

总结与体会
毕业论文终于结束了，我的大学生活也将随之结束。在这短短的一个多月里，我付出了许多，也收获了许多，回想起来还真有不少的体会。
首先要感谢我的导师彭昌荣老师，他在整个论文的完成中给了我极大的帮助和支持，他一丝不苟的治学态度，清晰的思路和认真、负责、勤快的工作作风深深的影响并感染了我，使我铭记于心并将永远的学习。在此，谨向恩师致以我最衷心的感谢。刚拿到题目时，我对课题充满了新鲜感，在导师的指导下我查阅了大量关于课题的文献，使我对课题的研究内容和国内外的研究进展有了更多的了解，也使我对自己的课题产生了浓厚的兴趣。
经过一个月的论文，我对自己也有了新的定位。那就是不能再把自己当做“小学生”了，干什么事情都要别人告诉我们该怎么做，我们是接受过高等教育的新时代的大学生，要完成时代赋予我们的使命，抓住机遇迎接挑战，就要锻炼自己发现问题和解决问题的能力，要养成主动学习的习惯，我认为这对我们以后的学习和工作是很重要的。
毕业论文是对我大学四年所学基础知识的专业知识的一次全面检验，通过撰写论文，我觉得对自己的语言组织能力，表达能力，沟通交际能力，运用所学知识的能力，分析问题并解决问题的能力都有所提高，也使我变的更加自信、成熟。
“团结、合作、谦虚”这三个词我的体会也比较深。做任何事包括做实验都不是孤立的，不是你“闭门造车”，而是一个需要和他人交往的过程。这就要求我们要团结，要有合作精神，要注意和他人的沟通，要谦虚，不懂就问所谓“知之为知之，不知为不知”。
总之，在整个论文完成的过程中，我体会到的是艰辛和收获的充实，感受到的是一种坚持不懈、契而不舍的科研精神。
最后，我再次向恩师和我的家人表示深深的感谢，我将继续努力，在以后的工作中，争取取得更好的成绩。
谢 辞
感谢彭昌荣老师对我的论文写作过程中的不足不厌其烦的细心指点。彭昌荣老师首先为我讲题，又对论文的格式的要求做了详细的讲解。当我迷茫于众多的资料时，他又为我提纲挈领，梳理脉络，使我确立本文的框架。论文的写作中，每周都得到了彭昌荣老师的指点。从框架的完善，到内容的扩充；从行文的用语，到格式的规范，彭昌荣老师都严格要求，力求完美。我再次为彭昌荣老师的付出表示由衷的感谢。

参考文献
[1] 胡英.物理化学(第四版)[M].北京:高等教育出版社,2002.
[2] 傅献彩.物理化学(第四版) [M].北京:高等教育出版社,1999.
[3] 许文林, 丁 平, 袁渭康. 硝基苯电解还原制备对氨基苯酚过程优化研究[J] . 化学工程, 1994, 22( 4) : 35~ 41.
[4] 刘大壮,徐海升,王安中,赵建宏. 连串反应工艺条件最优化J].化工高等教育,No2,34,1993.
[5] 徐海升王安中. 连串、平衡反应工程收率最优化[J]. 化学工程,Vol.17,No.2,18,1994.
[6] 金海林,藏福禄.应用最大值理化学反应的最佳温度分布[J].燕山石化,No.2,94,1989.
[7] 黄雪征, 张常群. 连续反应动力学的计算及过程的计算机模拟[J]. 计算机与应用化学，
Vol.21,No.3,211,2005.
[8] Gleich Walter. Continuous product ion of dialkanolamines[P] . DE3547328. 1994.
[9] 蒋旭峰, 曾崇余, 任晓乾等. 二异丙醇胺合成反应规律的研究[J].南京化工大学学报, 2000, 22( 2) : 48-511.
[10] 蒋旭峰. 高选择性合成二异丙醇胺反应动力学及新工艺研究[D]. 南京: 南京化工大学, 2000.
[11] 蒋旭峰, 曾崇余, 任晓乾, 黄永春. 二异丙醇胺合成反应动力学[J]. 南京:南京工业大学学报. Vol. 24,No. 1, Jan. 2002,24(1):102-104.
[12] 苏金明,阮沈勇,王永利编.MATLAB工程数学[M].北京:电子工业出版社,2005.
[13] Dadebo. S. A. and Mcauley. K. B. dynamic optimization of constrained chemical engineering problems using dynamic programming[j]. Computers chem. Engng Vol. 19 . No.5 . pp. 513- 525. 1995.
[14] Ray W. H., Advanced Process Control McGraw-Hill, New York (1981)..
[15] Renfro J. G., A. M. Morshedi and Osbjornsen. O. A. simultaneous optimization and solution of systems described by differential/algebraic equations. Computers chem. Engng 11,503-517 (1987).

附MATLAB程序
function PFRTempOpt_Bx
a0=3;%c10=3mol/dm3
b0=5; %c20=3mol/dm3
T=20:1:60;%温度℃
for i=1:length(T)
tend=10;%积分时间的最大值10min
[t,x,topt,x3max]=func_topt_x3max(a0,b0,tend,T(i));
tcal_topt(i)=topt;
x3max_cal(i)=x3max;
end
plot(T,x3max_cal,'-o')%打印各个温度T下的c3最大值
[S,E,G]=golden_point_opt_max(T(1),T(end),10^(-6),10^(-6),x3max_cal,T);
Topt=S(1)%c3达到最大值时的温度
x3max_global=S(2) %c3最大值
[t,x] = ode45(@CEquation,[0 10],[3 5 0]);%计算优化条件下：T=60℃，t=10min下的各个组分浓度分布，注意温度在Cequation中设定哈。
figure
plot(t,x(:,1),'-s')
hold on
plot(t,x(:,2),'-*')
hold on
plot(t,x(:,3),'-o')
function [t,x,topt,x3max]=func_topt_x3max(a0,b0,tend,T)
[t,x] = ode45(@(t,x) C1Equation(t,x,T),[0 tend],[a0 b0 0],T);
plot(t,x(:,1),'-o')
hold on
plot(t,x(:,2),'-s')
hold on
plot(t,x(:,3),'-d')
time=t;
xcal=x(:,3);
[S,E,G]=golden_point_opt_max(time(1),time(end),10^(-6),10^(-6),xcal,time);
topt=S(1);
x3max=S(2);
function dxdt = C1Equation(t,x,T) %带参数T的浓度方程
Ea1=9.141*10^4;
Ea2=9.130*10^4;
A1=7.967*10^12;
A2=4.682*10^12;
k1 = A1*exp(-Ea1/8.314/(T+273.15));
k2 = A2*exp(-Ea2/8.314/(T+273.15));
dxdt(1) =-k1*x(1)*x(2);
dxdt(2) =-k1*x(1)*x(2)-k2*x(2)*x(3);
dxdt(3) =k1*x(1)*x(2)-k2*x(2)*x(3);
dxdt = dxdt';
function dxdt = CEquation(t,x) %浓度方程
T=60;
Ea1=9.141*10^4;
Ea2=9.130*10^4;
A1=7.967*10^12;
A2=4.682*10^12;
k1 = A1*exp(-Ea1/8.314/(T+273.15));
k2 = A2*exp(-Ea2/8.314/(T+273.15));
dxdt(1) =-k1*x(1)*x(2);
dxdt(2) =-k1*x(1)*x(2)-k2*x(2)*x(3);
dxdt(3) =k1*x(1)*x(2)-k2*x(2)*x(3);
dxdt = dxdt';
function [S,E,G]=golden_point_opt_max(a,b,delta,epsilon,FXP,XP)%ÀëÉ¢µãÇó×î´óÖµ£ºÒÑÖªÏòÁ¿XP,ºÍFXP=f(XP)
%a,b-are the end poits of the interval
%delta-is the tolerance for the abscissas
%epsilon-is the tolerence for the ordinates
%Output-S=[p,yp]contains the abscissa p and the ordinates yp of the minimum
%E=[dp,dy] contains the error bounds for p and yp
%G-is an n¡Á4 matrix; kth row contains[ak,ck,dk,bk];the values of a,c,d and bat the kth iteration
r1=(sqrt(5)-1)/2;
r2=r1^2;
h=b-a;
[z]=interp1(XP,FXP,a,'spline');
ya=z;
[z]=interp1(XP,FXP,b,'spline');
yb=z;
c=a+r2*h;
d=a+r1*h;
[z]=interp1(XP,FXP,c,'spline');
yc=z;
[z]=interp1(XP,FXP,d,'spline');
yd=z;
k=1;
A(k)=a;B(k)=b;C(k)=c;D(k)=d;
while (abs(yb-ya)>epsilon)|(h>delta)
k=k+1;
if(yc>yd)
b=d;
yb=yd;
d=c;
yd=yc;
h=b-a;
c=a+r2*h;
[z]=interp1(XP,FXP,c,'spline');%f(c);
yc=z;
else
a=c;
ya=yc;
c=d;
yc=yd;
h=b-a;
d=a+r1*h;
[z]=interp1(XP,FXP,d,'spline');%f(d);
yd=z;
end
A(k)=a;B(k)=b;C(k)=c;D(k)=d;
end
dp=abs(b-a);
dy=abs(yb-ya);
p=a;
yp=ya;
if(yb>ya)
p=b;
yp=yb;
end
G=[A' C' D' B'];
S=[p,yp];
E=[dp dy];
附计算程序2
function CSTRTempOpt_lun%计算连串反应¦A-B-C使得B的浓度最大的温度分布
global T dt Ci0%温度，时间间隔，系统物质初始浓度
n=20%计算温度，时间间隔的数目
tspan =[0.00:1/n:1.00];
dt=1/n
Temp0 =398:-(398-298)/n:298
lb = ones(size(T))*298;
ub = ones(size(T))*398;
A = zeros(length(tspan),length(tspan));
%使得T1小于T2，T2小于T3…..Tn-1小于Tn，即使得温度由高%到低的%约束
b = zeros(length(tspan),1); %使得T1小于T2，T2小于T3…..Tn-1小于Tn，即使得温度由高%到低的约束
for i = 1:length(tspan)-1
A(i,i) = -1;
A(i,i+1) = 1;
end
A=A;
b=b;
Ci0=[1 0 0]
[T,fval,exitflag] = fmincon(@ObjFunc,Temp0,A,b,[],[],lb,ub)
plot(tspan,T,'-o')
xlabel('time/s')
ylabel('Temperature/¡æ')
stairs(tspan,T,'-')
TempCal=T%
x_dist = x_distCal(TempCal)%
figure
plot(tspan,x_dist(:,2),'-o')
xlabel('Time/s')
ylabel('c_2')
T=TempCal(1)%
[t,x] = ode45(@C1Equation,[0:dt/10:dt],Ci0)%
% ------------------------------------------------------------------
delta=0.0001;
epsilon=0.0000001;
for i=1:length(Temp0)
T=Temp0(i);
[t,x] = ode45(@C1Equation,[0:1/10:1],Ci0);
[S,E,G]=golden_point_opt_max(t(1),t(end),delta,epsilon,x(:,2),t);
time_opt(i)=S(1);
X2_max(i)=S(2);
end
X2_max=X2_max
[S,E,G]=golden_point_opt_max(Temp0(1),Temp0(end),delta,epsilon,X2_max,Temp0);
Temp_opt=S(1)
x2_max=S(2)
figure
plot(Temp0,X2_max,'-o')
hold on
plot(Temp_opt,x2_max,'-rs')%c2最大时的温度
xlabel('Temperature/K')
ylabel('c_2_m_a_x')
grid on

%------------------------------------------------------------------
function f = ObjFunc(Temp)% 目标函数
global T dt Ci0
x0=Ci0;
for i=1:length(Temp)
T=Temp(i);
[t,x] = ode45(@C1Equation,[0:dt/10:dt],x0);
x0=x(end,:);
end
f=-x(end,2); % 目标函数-使得B的浓度最大
%----------------------------
function x_dist = x_distCal(TempCal)%
global T dt Ci0
x0=Ci0;
for i=1:length(TempCal)
T=TempCal(i);
[t,x] = ode45(@C1Equation,[0:dt/10:dt],x0);
x0=x(end,:);
x_dist(i,:)=x(end,:);
end
function dxdt = C1Equation(t,x) %微分方程
global T
Ea1=2500;
Ea2=5000;
A1=4*10^3;
A2=620*10^3;
k1=A1*exp(-Ea1./T);
k2=A2*exp(-Ea2./T);
dxdt(1) =-k1*x(1)*x(1);
dxdt(2) =k1*x(1)*x(1)-k2*x(2);
dxdt(3) =k2*x(2);
dxdt = dxdt';
26