华南理工大学化工学院2001级博士生课程设计

PAGE华南理工大学化工学院2001级博士生课程设计超临界CO2萃取紫苏油过程的计算机模拟设计者:文震指导老师:钱宇2001年12月目录摘要………………………1Abstract…………………2一设计背景……………3二设计原理……………3BP神经网络原理……………………3BP神经网络算法…………………4BP算法程序流程图…………………6网络参数的选择……………………6三实验部分…………61实验原料……………62实验试剂与设备……………………63实验方法与流程………7四结果与讨论………71神经网络的训练……………………82网络有效性的验证…………………93网络预测结果………10结论………………………11参考文献…………………12附录…………13PAGE1-图1BP神经网络拓扑结构BPN的学习由以下四个过程组成：学习模式输入神经网络后，神经元的激活值由输入层经隐含层，在输出层各神经元获得响应的“模式顺传播”过程；然后按减小希望输出与实际输出的误差方向，从输出层向隐含层再向输入层逐层修正各连接权值的“误差逆传播”过程；由“模式顺传播”过程与“误差逆传播”过程的交替进行的网络学习训练过程；网络全局误差趋向极小的学习收敛过程。当全局误差小于设定值或学习次数大于设定值时，均应停止训练。然后对训练好的网络进行性能测试，即转入网络预测阶段[3,4]。2误差反向传播（BP）算法BP学习算法适用于层次型神经网络。图1即表示一个三层学习网络，其中包括输入层、输出层和两者之间的隐含层。同层神经元间没有连接，相邻层间有全连接，由低层指向高层。设层次网络共有m+1层，输入层为第0层，输出层为第m层。令Weq\o(\s\up6(ｎ),\s\do6(ｊｋ))表示第n+1层中第j神经元与第n中第k单元间的连接权，则神经元输出的计算公式为：Xeq\o(\s\up6(ｎ＋１),\s\do6(ｊ))=f(Weq\o(\s\up6(ｎ),\s\do6(ｊｋ))Xeq\o(\s\up6(ｎ),\s\do6(ｋ)))n=0,1,…,m(1)式中f表示S型函数。本文选择的S型函数为：f(u)=1/1(1+exp(-u))当选取定L个训练样本后，调整Weq\o(\s\up6(ｎ),\s\do6(ｊｋ))的目的是使输出层神经元的实际输出和应有的已知输出的平方误差尽可能小，即使得E=(deq\o(\s\up6(ｐ),\s\do6(ｊ))-Xeq\o(\s\up6(ｐｍ),\s\do6(ｊ)))2(2)最小，式中：deq\o(\s\up6(ｐ),\s\do6(ｊ))为第p个样本输入时网络中第j个单元的已知输出值;Xeq\o(\s\up6(ｐｍ),\s\do6(ｊ))为第p个样本输入时，第m层第j个单元的实际输出值。调整Weq\o(\s\up6(ｎ),\s\do6(ｊｋ))的公式为：Weq\o(\s\up6(ｐｎ),\s\do6(ｊｋ))=Weq\o(\s\up6(ｐ－１ｎ),\s\do6(ｋ))+△Weq\o(\s\up6(ｐｎ),\s\do6(ｊｋ))(3)为了减小E，采用使E沿梯度方向下降的方式，即使得修正量△Weq\o(\s\up6(ｐｎ),\s\do6(ｊｋ))与-δ(E)/δ(Weq\o(\s\up6(ｐｎ),\s\do6(ｊｋ)))成正比。由式（1）)、（2）知，对于p=1,2,…L，可写出一般式(-δ(E)/δ(Weq\o(\s\up6(ｐｎ),\s\do6(ｊｋ))))*2（dj-Xeq\o(\s\up6(ｎ＋１),\s\do6(ｊ)))(δXeq\o(\s\up6(ｎ＋１),\s\do6(ｊ))/δ(Weq\o(\s\up6(ｐｎ),\s\do6(ｊｋ))))=2（dj-Xeq\o(\s\up6(ｎ＋１),\s\do6(ｊ)))(δf/δ)*(δ/δ(Weq\o(\s\up6(ｐｎ),\s\do6(ｊｋ))))=2（dj-Xeq\o(\s\up6(ｎ＋１),\s\do6(ｊ)))*Xeq\o(\s\up6(ｎ＋１),\s\do6(ｊ))(1-Xeq\o(\s\up6(ｎ＋１),\s\do6(ｊ)))*Xeq\o(\s\up6(ｎ),\s\do6(ｋ))当n+1=m时，令δeq\o(\s\up6(ｍ),\s\do6(ｊ))=(dj-Xeq\o(\s\up6(ｍ),\s\do6(ｊ)))Xeq\o(\s\up6(ｍ),\s\do6(ｊ))(1-Xeq\o(\s\up6(ｍ),\s\do6(ｊ)))(4)由于dj已知，于是δeq\o(\s\up6(ｍ),\s\do6(ｊ))可以求出。因△Weq\o(\s\up6(ｍ－１),\s\do6(ｊｋ))=ηδeq\o(\s\up6(ｍ),\s\do6(ｊ))Xeq\o(\s\up6(ｍ－１),\s\do6(ｋ)),0≤η≤1(5)此处δeq\o(\s\up6(ｍ),\s\do6(ｊ))是输出层的误差。以后各层误差是将上一层误差通过网络往下传播而得，即δeq\o(\s\up6(ｎ),\s\do6(ｊ))=Xeq\o(\s\up6(ｎ),\s\do6(ｊ))(1-Xeq\o(\s\up6(ｎ),\s\do6(ｊ)))δeq\o(\s\up6(ｎ＋１),\s\do6(ｋ))Weq\o(\s\up6(ｎ),\s\do6(ｋｊ))(6)相应的Wjk调整值为△Weq\o(\s\up6(ｎ),\s\do6(ｊｋ))=ηδeq\o(\s\up6(ｎ＋１),\s\do6(ｊ))Xeq\o(\s\up6(ｎ),\s\do6(ｋ))(7)因此，BP学习算法的具体步骤是：1、赋予各Wjk以较小的随机非零初始值；2、输入一个样本和它已知的输出值；3、按公式(1)计算出实际的输出值；4、按公式(4)至(7)修改连接权值；5、转入第2步，循环运用L个样本作为输入，直至Weq\o(\s\up6(ｎ),\s\do6(ｊｋ))均稳定为止。3BP算法程序流程图4网络参数的选定本文采用的BPN由输入层，一层隐含层，输出层组成。考虑到压力与温度是影响超临界CO2萃取紫苏油过程的主要因素，输入变量为两个，分别对应压力与温度，输出变量为1个，即此温度压力下所对应的溶解度。隐含层的节点数根据全局误差最小原则确定。二实验部分1实验原料紫苏子：产自贵阳花溪，干燥后粉碎至40目2实验试剂与设备CO2：贵阳都拉综合厂，食品级，>99.5%超临界CO2萃取装置：中科院地化所超临界中心研制，萃取器：4L1个，分离器：2L1个，系统最高使用压力32Mpa，装置如图2。CO2钢瓶2-过滤器3-冷机4-高压计量泵5-换热器6-萃取器7、8-分离器1、29-转子流量计10-累计流量计图2超临界CO2萃取装置流程图3实验方法与流程称取一定量紫苏子放入萃取罐中，CO2从钢瓶出来后，经冷机冷凝，由高压计量泵加压至设定压力，由底部进入萃取器。循环水加热预热器、萃取器和分离器至所需温度，调节CO2流量进行萃取。SC-CO2携带着紫苏油进入分离器，升温减压后，紫苏油从SC-CO2中析出，由分离器下部阀门放出，用接受瓶接受。在实验中，要保持尽可能小的CO2流速，使紫苏油在SC-CO2中溶解状态接近平衡饱和状态，这种因设备限制而采取近似的方法来测定待萃物的溶解度的方法被普遍使用[5]。紫苏油在SC-CO2中溶解度S按下式计算：S＝M/Δt×FM－Δt时间内从分离器接受紫苏油的质量（g），F－CO2流速(Kg/h)三结果与讨论紫苏油在SC-CO2中溶解度测定结果如表1。由表1可见，溶解度是温度与压力的函数，温度升高，SC-CO2密度下降，溶解度降低，压力增大，SC-CO2密度增加，溶解度增大。表1紫苏油在SC-CO2中溶解度（g/Kg）温度（℃）萃取压力（Mpa）151821242730350.791.753.02*4.575.736.47450.411.18*2.564.185.496.34550.280.882.103.645.196.08*650.19*0.631.653.254.845.97750.080.431.252.73*4.725.701紫苏油在SC-CO2中溶解度的神经网络的训练1．1隐含层节点数的确定目前隐含层节点个数的选择尚无理论规则，本文利用误差最小原则来确定节点数。网络误差与隐节点数的关系如图3。图3网络误差与隐节点数的关系由图3可以看出，在节点数≤5时，网络误差随节点数增加而下降，当节点数＞5时，网络误差变化不大，节点为9时，误差反而略为增加，而且节点数越多，网络计算量越大，计算速度越慢，因此，隐含层节点数选择5为宜。1．2迭代次数的确定我们将试验样本随机地分为训练集与测试集，其中训练集25组，测试集5组，测试集如表1打“*”数据。网络误差与迭代次数的关系如图4。图4网络误差与迭代次数的关系图4显示：在迭代次数＜2000时，训练集与测试集的网络误差急剧下降，但随着训练的继续进行，迭代次数增加，网络误差的下降趋势变缓，当迭代次数＞4000后，测试集的误差反而略有增加。这可能是因为出现过拟合的缘故，过拟合产生是由于网络模型刻意去契合个别样本所致，此时，网络对参加学习的样本可以给出较好的拟合，但对于不参加训练的样本，拟合偏差反而增大，因此不管此时训练集的误差是否继续下降，均应停止迭代。综合训练集与测试集的网络误差考虑，迭代次数应定为5000为宜。2网络有效性的验证对于测试集，网络输出与理想输出（实际值）的对比如图5。图5网络输出与理想输出的对比图5显示：网络的输出值与实际值非常接近，最大相对误差为9.30％，最小相对误差为0.04％，整个测试集的误差仅为3.27％，这说明网络具有较好的模拟效果。此外，此时训练集的平均误差为5.40％，比测试集稍大，这是因为训练集中含有因试验误差所致的个别异常点，使得网络输出与这些异常值的误差较大，从而增加了整个训练集的平均相对偏差。这也同时也表明，神经网络具有一定的泛化能力，个别异常样本出现不影响整个网络的预报精度，说明紫苏油溶解度的神经网络模型是有效的。3网络预测结果利用已建立好的网络模型，我们对萃取参数为压力15-30Mpa，步长为1Mpa，温度为35-75℃下的溶解度进行了预测。预测结果如图6。图6紫苏油在SC-CO2中的溶解度规律图6清楚地显示了紫苏油在SC-CO2中的溶解规律：压力是主要因素，温度是次要因素。在压力较低阶段，紫苏油几乎不溶于SC-CO2中，在压力为22－25MPa阶段，溶解度的增幅最大，当压力越过25MPa后，这种增幅明显趋缓。温度对溶解度则呈负影响，即温度增加，溶解度下降。在较低压力下，由于溶解度的绝对值较小，温度的影响也较小。在压力影响最显著阶段，温度的影响也最大，高温与低温之间的溶解度差值最大。随着压力的继续上升，这种差值变小，说明此时压力已开始接近转折压力，整个萃取过程将转为温度正效应阶段。受设备所限，无法达到转折压力（据文献报道，油脂类物质的转折压力在40Mpa左右[6]），但这也表明神经网络的模拟结果是符合溶解度变化的实际规律的。四结论利用间歇萃取装置，采取静态法测定并计算了紫苏油在SC-CO2中的溶解度。建立并验证了紫苏油在SC-CO2中溶解度的神经网络模型，该模型能较好地拟合实际溶解度。（3）利用已建立的溶解度模型，并不需要了解物质在SC-CO2中的溶解机理，即可预测未知条件下的溶解度，反映溶解规律，从而可以为萃取工艺的设计与优化，以及SCFE装置的工业放大数据的计算提供理论依据。参考文献[1]Liong.K.K.,FosterN.R.andTingS.S.T.,Solubilityoffattyacidestersinsupercriticalcarbondioxide,Ind.Eng.Chem.Res.,1992,31(1):400-404[2]BambergerT.,EricksonJ.C.,CooneyC.L.,etal.,Measurementandmodelpredictionofsolubilitiesofpurefattyacids,puretriglyceridesandmixturesoftriglyceridesinsupercriticalcarbondioxide,J.Chem.Eng.Data.,1988,33(3):327-333.[3]陈念贻，钦佩，陈瑞亮，模式识别方法在化学化工中的应用，科学出版社，2000，北京，71-101[4]阎平凡，张长水，人工神经网络与模拟进化计算，清华大学出版社，2000，北京，10-37[5]马海乐，吴守一，陈均等，小麦胚芽油在超临界CO2中溶解度的试验研究，中国粮油学报，1997，12（1）：38-47[6]NilssonW.B.,Ganglitz.Jr.E.J.andHudsonJ.K.,SolubilitiesofmethylOleate,Oleicacid,Oleylglycerols,andOleylglycerolmixturesinsupercriticalcarbondioxide,J.Am.Oilchem.Soc.,1991,68(2):87-91.附录程序源代码一、模块源代码1、定义变量OptionExplicit'定义输入输出的文件名Publicfile_inputAsStringPublicfile_weightAsString'定义输入和输出的数据个数Publicinput_numberAsIntegerPublicoutput_numberAsInteger'定义每层的结点数Publicnodes_of_layer0%,nodes_of_layer1%,nodes_of_layer2%Publiclearning_rateAsSingle'输入数值和输出数值;'原始数据Publicinputdata()AsSinglePublicoutputdata()AsSingle'归一化以后的数据Publicnetinput()AsSinglePublicnetoutput()AsSingle'各层权重Publicweight_10()AsSinglePublicweight_21()AsSingle'各层的输出值Publiclayer0_out()AsSinglePubliclayer1_out()AsSinglePubliclayer2_out()AsSinglePublicerror1()AsSinglePublicerror2()AsSinglePublicprogram_stateAsStringPublicregulate_inputAsSinglePublicregulate_outputAsSingle2、重新设置'本子程序的目的:重置所有的数组,重新开始计算Subrefresh_net()Dimi%,j%ReDiminputdata(nodes_of_layer0)AsSingleReDimoutputdata(nodes_of_layer0)AsSingleReDimnetinput(nodes_of_layer0)AsSingleReDimnetoutput(nodes_of_layer0)AsSingleReDimweight_10(nodes_of_layer1,nodes_of_layer0)AsSingleReDimweight_21(nodes_of_layer2,nodes_of_layer1)AsSingleReDimlayer0_out(nodes_of_layer0)AsSingleReDimlayer1_out(nodes_of_layer1)AsSingleReDimlayer2_out(nodes_of_layer2)AsSingleReDimerror1(nodes_of_layer1)AsSingleReDimerror2(nodes_of_layer2)AsSingle3、权重随机赋值'给各层的权重赋随机数RandomizeFori=1Tonodes_of_layer1Forj=1Tonodes_of_layer0weight_10(i,j)=RndNextjNextiFori=1Tonodes_of_layer2Forj=1Tonodes_of_layer1weight_21(i,j)=RndNextjNextiEndSub4、计算正向输出结果'该子程序用于正向计算输出结果Subcaculate()Dimi%,j%DimftempAsSingle'caculatelayer0Fori=1Tonodes_of_layer0layer0_out(i)=netinput(i)Nexti'caculatelayer1Fori=1Tonodes_of_layer1ftemp=0Forj=1Tonodes_of_layer0ftemp=ftemp+weight_10(i,j)*layer0_out(j)Nextjftemp=1/(1+Exp(-ftemp))layer1_out(i)=ftempNexti'caculatelayer2Fori=1Tonodes_of_layer2ftemp=0Forj=1Tonodes_of_layer1ftemp=ftemp+weight_21(i,j)*layer1_out(j)Nextjftemp=1/(1+Exp(-ftemp))layer2_out(i)=ftempNextiEndSub5、调整网络权重'本子程序的作用是:对网络的权重进行调整Subadjust()Dimi%,j%DimftempAsSingle'调整输出层的权重Fori=1Tonodes_of_layer2error2(i)=(netoutput(i)-layer2_out(i))*layer2_out(i)*(1-layer2_out(i))NextiFori=1Tonodes_of_layer2Forj=1Tonodes_of_layer1weight_21(i,j)=weight_21(i,j)+learning_rate*error2(i)*layer1_out(j)NextjNexti'调整隐含层的权重Fori=1Tonodes_of_layer1ftemp=0Forj=1Tonodes_of_layer2ftemp=ftemp+weight_21(j,i)*error2(j)Nextjerror1(i)=ftemp*layer1_out(i)*(1-layer1_out(i))NextiFori=1Tonodes_of_layer1Forj=1Tonodes_of_layer0weight_10(i,j)=weight_10(i,j)+learning_rate*layer0_out(j)*error1(i)NextjNextiEndSub6、数据的保存Subsave_net()Dimi%,j%OnErrorGoTo1000Openfile_weightForOutputAs#1Print#1,input_numberPrint#1,output_numberPrint#1,nodes_of_layer0,nodes_of_layer1,nodes_of_layer2Print#1,learning_ratePrint#1,regulate_inputPrint#1,regulate_outputFori=1Tonodes_of_layer1Forj=1Tonodes_of_layer0Print#1,weight_10(i,j);Nextj,iPrint#1,Fori=1Tonodes_of_layer2Forj=1Tonodes_of_layer1Print#1,weight_21(i,j);Nextj,iClose#1ExitSub1000:Close#1MsgBox"无法写入配置文件,无法记录权值",vbOKOnly,""Err.ClearEndSub7、自定义函数Functionload_net()AsStringDimi%,j%OnErrorGoTo1000Openfile_weightForInputAs#1Input#1,input_numberInput#1,output_numberInput#1,nodes_of_layer0,nodes_of_layer1,nodes_of_layer2Input#1,learning_rateInput#1,regulate_inputInput#1,regulate_output'重新定义数组refresh_netFori=1Tonodes_of_layer1Forj=1Tonodes_of_layer0Input#1,weight_10(i,j)Nextj,iFori=1Tonodes_of_layer2Forj=1Tonodes_of_layer1Input#1,weight_21(i,j)Nextj,iClose#1load_net="OK"ExitFunction1000:Close#1MsgBox"无法打开配置文件,重新开始学习",vbOKOnly,""Err.ClearEndFunction二、网络初始化及参数设置源代码OptionExplicitPrivateSubForm_Load()Text1.Text=nodes_of_layer1Text2.Text=Format(learning_rate,"0.00")UpDown1.Value=Text1.TextText3.Text=input_numberText4.Text=output_numberUpDown2.Value=Text3.TextUpDown3.Value=Text4.TextText5.Text=Form1.learn_timesText6.Text=Form1.max_error*100IfForm1.method=0ThenOption1.Value=TrueElseOption2.Value=TrueEndSubPrivateSubImage2_Click()UnloadMeEndSubPrivateSubImage5_Click()IfText2.Text<=0OrText2.Text>1ThenMsgBox"学习效率的数值只能在0到1之间，请改正",vbOKOnly+vbExclamation,"提示栏！"Text2.SetFocusEndIfForm1.learn_times=Text5.TextForm1.max_error=Text6.Text/100IfOption1.Value=TrueThenForm1.method=0ElseForm1.method=1learning_rate=Text2.TextIfText1.Text<>nodes_of_layer1OrText3.Text<>input_numberOr_Text4.Text<>output_numberThen'Orlearning_rate<>Text2.TextDimresponse%response=MsgBox("网络参数已被改变,是否要重置网络?选定后按[返回]键返回。",vbYesNo+vbExclamation,"提示栏！")Ifresponse=vbYesThenForm1.total_times=0input_number=Text3.Textoutput_number=Text4.Textnodes_of_layer1=Text1.Text'考虑到阀值的情况IfCheck1.Value=0Thennodes_of_layer0=Text3.Textnodes_of_layer2=Text4.TextEndIfrefresh_netIfForm1.data_okThensave_netEndIfEndIfEndSubPrivateSubImage6_Click()IfForm1.data_ok=TrueThenrefresh_netsave_netEndIfEndSubPrivateSubText5_Change()Option1.Value=TrueEndSubPrivateSubText6_Change()Option2.Value=TrueEndSubPrivateSubUpDown1_Change()Text1.Text=UpDown1.ValueEndSubPrivateSubUpDown2_Change()Text3.Text=UpDown2.ValueEndSubPrivateSubUpDown3_Change()Text4.Text=UpDown3.ValueEndSub三、学习训练窗口源代码OptionExplicitDimarray1()AsSingleDimarray2()AsSingleDimgroups%Publictotal_timesAsLongPublicdata_okAsBooleanPubliclearn_timesAsLongPublicmax_errorAsSinglePublicmethod%PrivateSubForm_Load()nodes_of_layer0=26nodes_of_layer1=8nodes_of_layer2=1input_number=1output_number=1max_error=0.01learn_times=10000method=0learning_rate=0.9EndSubFunctionload_data(file_dAsString,array1()AsSingle,array2()AsSingle)AsIntegerDimi%,j%OnErrorGoTo1000Openfile_dForInputAs#1Dimgroups%,totalAsLongReDimarray1(1000,input_number)ReDimarray2(1000,output_number)groups=1total=0regulate_input=0regulate_output=0DoWhileNotEOF(1)Fori=1Toinput_numberInput#1,array1(groups,i)Ifregulate_input<Abs(array1(groups,i))Thenregulate_input=Abs(array1(groups,i))total=total+1NextiForj=1Tooutput_numberInput#1,array2(groups,j)Ifregulate_output<Abs(array2(groups,j))Thenregulate_output=Abs(array2(groups,j))total=total+1Nextjgroups=groups+1Loop1000:IfEOF(1)Thenregulate_input=regulate_input*2regulate_output=regulate_output*2groups=total/(input_number+output_number)load_data=groupsClose#1ExitFunctionElseClose#1MsgBox"无法打开数据文件",vbOKOnly,""Err.ClearEndIfEndFunctionPrivateSubImage1_Click()Ifdata_okThensave_netEndSubPrivateSubImage2_Click()Dimi%,j%Dimx%IfNotdata_okThenMsgBox"请先打开数据文件!",vbOKOnly,""ExitSubEndIfDimftempAsSingleSelectCasemethodCase0:Forx=1Tolearn_timesftemp=0Fori=1TogroupsForj=1Toinput_numbernetinput(j)=array1(i,j)/regulate_inputNextjForj=1Tooutput_numbernetoutput(j)=array2(i,j)/regulate_outputNextjcaculateForj=1Tooutput_numberftemp=ftemp+Abs((layer2_out(j)-netoutput(j))/netoutput(j))NextjadjustNextiDoEventsImage5.Visible=Truetotal_times=total_times+1Text2.Text=total_timesftemp=ftemp/(output_number*groups)Text1.Text=Format(ftemp*100,"0.000000000")&"%"NextxCase1:Doftemp=0Fori=1TogroupsForj=1Toinput_numbernetinput(j)=array1(i,j)/regulate_inputNextjForj=1Tooutput_numbernetoutput(j)=array2(i,j)/regulate_outputNextjcaculateForj=1Tooutput_numberftemp=ftemp+Abs((layer2_out(j)-netoutput(j))/netoutput(j))NextjadjustNextiftemp=ftemp/(output_number*groups)Text2.Text=total_timesText1.Text=Format(ftemp*100,"0.000000")&"%"I