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

1、华南理工大学化工学院2001级博士生课程设计超临界CO2萃取紫苏油过程的计算机模拟设计者: 文 震指导老师: 钱 宇2001年12月目 录摘要1Abstract2一 设计背景3二 设计原理31 BP神经网络原理32 BP神经网络算法 43 BP算法程序流程图64 网络参数的选择 6三 实验部分 61 实验原料 62 实验试剂与设备 63实验方法与流程 7四 结果与讨论 71 神经网络的训练 82 网络有效性的验证 93 网络预测结果 10结论 11参考文献12附录13超临界CO2萃取紫苏油过程的计算机模拟摘 要本文测定了紫苏油在超临界CO2（SC-CO2）中的溶解度，并利用误差逆传播（BP）神

2、经网络对溶解度数据进行了拟合。通过对萃取参数与溶解度的关系进行训练学习，实现网络结构的优化，建立了紫苏油在SC-CO2溶解度的网络模型，并将该网络用于未知萃取参数下溶解度的预测，得到了较高的预测精度。结果表明，该方法可作为预测物质在SC-CO2溶解度的一种有效手段。关键词： 溶解度 紫苏油 神经网络 模型 超临界CO2Simulation of Extraction Process of Perilla Seed Oil with Supercritical Carbon DioxideAbstractIn this article, the solubility of Perilla see

3、d oil in supercritical carbon dioxide (SC-CO2) was measured, then the solubility was simulated by Back Propagation (BP) neural network. By means of learning and training with relation of extraction parameter and solubility , the neural model of the solubility of Perilla Seed Oil in supercritical car

4、bon dioxide was established, then solubility under unknown extraction parameter was predicted with this model, the precision of the prediction is fairly high. The results showed this method can be used as an effective way to predict the solubility of matter in SC-CO2.Key words: solubility Perilla se

5、ed oil neural network model supercritical carbon dioxide 一 设计背景了解和预测待萃物在超临界CO2流体中的溶解度及其随温度、压力等影响因素的变化规律, 是进行超临界CO2萃取工艺设计和设备计算的基础。而目前待萃物的溶解度主要由实验所测，受主观影响因素较大。近年来，研究者尝试用建立数学模型的方法来计算纯物质的溶解度，而在超临界状态下物质的理化参数极其匮乏，同时由于待萃物大多数为组成复杂的天然产物，很难利用现成的纯物质的溶解度公式建立数学模型，因此寻求一种能预测混合物溶解度方法是很有必要的 1, 2 。人工神经网络(Artificial N

6、eural Network，简称ANN)是80年代中后期迅速发展起来的研究领域, 主要用于函数逼近与分类识别优化，特别适用于复杂系统的数学模型的建立 3,4 。本文选取超临界CO2萃取紫苏油这一较为典型的动态、高度非线性的过程，拟利用误差逆传播(Back Propagation，BP)神经网络对紫苏油在SC-CO2中溶解度进行拟合与建模，来预测未知条件下的溶解度。期望能通过利用试验结合数学模型的方法，为超临界CO2萃取紫苏油工艺的设计与优化提供理论依据。二 设计原理1 BP神经网络原理BP神经网络（BPN）是一种具有三层神经元的阶层网络，不同阶层的神经元之间实现权重连接，其拓扑结构如图1。图中

7、，m、p、q分别是输入层、隐含层、输出层，Wi j为输入层与隐含层的连接权值，Vj t为隐含层至输出层的连接权值。Ak输入模式向量，Yk为输出模式向量。隐含层j=1,2, ,p输入层Ak=(a1k,a2k.,amk)I=1,2, , mWi,j输出层Yk=(y1k,y2k.yqk)t=1,2, , qVj,t图1 BP神经网络拓扑结构BPN的学习由以下四个过程组成：学习模式输入神经网络后，神经元的激活值由输入层经隐含层，在输出层各神经元获得响应的“模式顺传播”过程；然后按减小希望输出与实际输出的误差方向，从输出层向隐含层再向输入层逐层修正各连接权值的“误差逆传播”过程；由“模式顺传播”过程与“

8、误差逆传播”过程的交替进行的网络学习训练过程；网络全局误差趋向极小的学习收敛过程。当全局误差小于设定值或学习次数大于设定值时，均应停止训练。然后对训练好的网络进行性能测试，即转入网络预测阶段 3,4 。2 误差反向传播（BP）算法 BP学习算法适用于层次型神经网络。图1即表示一个三层学习网络，其中包括输入层、输出层和两者之间的隐含层。同层神经元间没有连接，相邻层间有全连接，由低层指向高层。 设层次网络共有m+1层，输入层为第0层，输出层为第m层。令W表示第n+1层中第j神经元与第n中第k单元间的连接权，则神经元输出的计算公式为：X=f(WX) n=0,1,m (1)式中f表示S型函数。本文选择

9、的S型函数为：f(u)=1/1(1+exp(-u) 当选取定L个训练样本后，调整W的目的是使输出层神经元的实际输出和应有的已知输出的平方误差尽可能小，即使得E=(d-X) 2 (2)最小，式中： d为第p个样本输入时网络中第j个单元的已知输出值; X为第p个样本输入时，第m层第j个单元的实际输出值。调整W的公式为： W=W +W (3) 为了减小E，采用使E沿梯度方向下降的方式，即使得修正量W与-(E)/(W)成正比。由式（1）)、（2）知，对于p=1,2,L，可写出一般式 (-(E)/(W) *2（dj-X)(X/(W)= 2（dj-X) (f/)*(/(W)= 2（dj-X)*X(1-X)

10、*X当n+1=m时，令 =(dj-X)X(1-X) (4)由于dj已知，于是可以求出。因 W=X , 0 1 (5)此处是输出层的误差。 以后各层误差是将上一层误差通过网络往下传播而得，即 =X(1-X)W (6)相应的Wjk调整值为 W= X (7) 因此，BP学习算法的具体步骤是：1、赋予各Wjk以较小的随机非零初始值；2、输入一个样本和它已知的输出值；3、按公式(1)计算出实际的输出值；4、按公式(4)至(7)修改连接权值；5、转入第2步，循环运用L个样本作为输入，直至W均稳定为止。3 BP算法程序流程图4 网络参数的选定本文采用的BPN由输入层，一层隐含层，输出层组成。考虑到压力与温度

11、是影响超临界CO2萃取紫苏油过程的主要因素，输入变量为两个，分别对应压力与温度，输出变量为1个，即此温度压力下所对应的溶解度。隐含层的节点数根据全局误差最小原则确定。二 实验部分1 实验原料紫苏子：产自贵阳花溪，干燥后粉碎至40目2 实验试剂与设备CO2：贵阳都拉综合厂，食品级，99.5%超临界CO2萃取装置：中科院地化所超临界中心研制，萃取器：4L1个，分离器：2L 1个，系统最高使用压力32Mpa，装置如图2。1- CO2钢瓶 2-过滤器 3-冷机 4-高压计量泵 5-换热器 6-萃取器 7、8-分离器1、2 9-转子流量计 10- 累计流量计图2 超临界CO2萃取装置流程图3 实验方法与

12、流程称取一定量紫苏子放入萃取罐中，CO2从钢瓶出来后，经冷机冷凝，由高压计量泵加压至设定压力，由底部进入萃取器。循环水加热预热器、萃取器和分离器至所需温度，调节CO2流量进行萃取。SC-CO2携带着紫苏油进入分离器，升温减压后，紫苏油从SC-CO2中析出，由分离器下部阀门放出，用接受瓶接受。在实验中，要保持尽可能小的CO2流速，使紫苏油在SC-CO2中溶解状态接近平衡饱和状态，这种因设备限制而采取近似的方法来测定待萃物的溶解度的方法被普遍使用 5 。紫苏油在SC-CO2中溶解度S按下式计算：SM / tFMt时间内从分离器接受紫苏油的质量（g），FCO2流速(Kg / h)三 结果与讨论紫苏油

13、在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中溶解度的神经网络的训练11 隐含层节点数

14、的确定目前隐含层节点个数的选择尚无理论规则，本文利用误差最小原则来确定节点数。网络误差与隐节点数的关系如图3。图3 网络误差与隐节点数的关系由图3可以看出，在节点数5时，网络误差随节点数增加而下降，当节点数5时，网络误差变化不大，节点为9时，误差反而略为增加，而且节点数越多，网络计算量越大，计算速度越慢，因此，隐含层节点数选择5为宜。12 迭代次数的确定我们将试验样本随机地分为训练集与测试集，其中训练集25组，测试集5组，测试集如表1打“*”数据。网络误差与迭代次数的关系如图4。图4 网络误差与迭代次数的关系图4显示：在迭代次数2000时，训练集与测试集的网络误差急剧下降，但随着训练的继续进行

15、，迭代次数增加，网络误差的下降趋势变缓，当迭代次数4000后，测试集的误差反而略有增加。这可能是因为出现过拟合的缘故，过拟合产生是由于网络模型刻意去契合个别样本所致，此时，网络对参加学习的样本可以给出较好的拟合，但对于不参加训练的样本，拟合偏差反而增大，因此不管此时训练集的误差是否继续下降，均应停止迭代。综合训练集与测试集的网络误差考虑，迭代次数应定为5000为宜。2 网络有效性的验证对于测试集，网络输出与理想输出（实际值）的对比如图5。图5 网络输出与理想输出的对比图5显示：网络的输出值与实际值非常接近，最大相对误差为9.30，最小相对误差为0.04，整个测试集的误差仅为3.27，这说明网络

16、具有较好的模拟效果。此外，此时训练集的平均误差为5.40，比测试集稍大，这是因为训练集中含有因试验误差所致的个别异常点，使得网络输出与这些异常值的误差较大，从而增加了整个训练集的平均相对偏差。这也同时也表明，神经网络具有一定的泛化能力，个别异常样本出现不影响整个网络的预报精度， 说明紫苏油溶解度的神经网络模型是有效的。3 网络预测结果利用已建立好的网络模型，我们对萃取参数为压力15-30Mpa，步长为1Mpa，温度为35-75下的溶解度进行了预测。预测结果如图6。图6 紫苏油在SC-CO2中的溶解度规律图6清楚地显示了紫苏油在SC-CO2中的溶解规律：压力是主要因素，温度是次要因素。在压力较低

17、阶段，紫苏油几乎不溶于SC-CO2中，在压力为2225MPa阶段，溶解度的增幅最大，当压力越过25MPa后，这种增幅明显趋缓。温度对溶解度则呈负影响，即温度增加，溶解度下降。在较低压力下，由于溶解度的绝对值较小，温度的影响也较小。在压力影响最显著阶段，温度的影响也最大，高温与低温之间的溶解度差值最大。随着压力的继续上升，这种差值变小，说明此时压力已开始接近转折压力，整个萃取过程将转为温度正效应阶段。受设备所限，无法达到转折压力（据文献报道，油脂类物质的转折压力在40Mpa左右6），但这也表明神经网络的模拟结果是符合溶解度变化的实际规律的。四 结论（1） 利用间歇萃取装置，采取静态法测定并计算了

18、紫苏油在SC-CO2中的溶解度。（2） 建立并验证了紫苏油在SC-CO2中溶解度的神经网络模型，该模型能较好地拟合实际溶解度。（3） 利用已建立的溶解度模型，并不需要了解物质在SC-CO2中的溶解机理，即可预测未知条件下的溶解度，反映溶解规律，从而可以为萃取工艺的设计与优化，以及SCFE装置的工业放大数据的计算提供理论依据。参考文献1 Liong. K .K., Foster N. R. and Ting S.S.T., Solubility of fatty acid esters in supercritical carbon dioxide, Ind. Eng. Chem. Res.,

19、1992, 31(1): 400-404 2 Bamberger T., Erickson J.C., Cooney C.L., et al., Measurement and model prediction of solubilities of pure fatty acids, pure triglycerides and mixtures of triglycerides in supercritical carbon dioxide, J. Chem. Eng. Data., 1988, 33(3): 327-333.3 陈念贻，钦佩，陈瑞亮，模式识别方法在化学化工中的应用，科学出版

20、社，2000，北京，71-1014 阎平凡，张长水，人工神经网络与模拟进化计算，清华大学出版社，2000，北京，10-375马海乐，吴守一，陈均等，小麦胚芽油在超临界CO2中溶解度的试验研究，中国粮油学报，1997，12（1）：38-476 Nilsson W. B., Ganglitz. Jr. E. J. and Hudson J.K., Solubilities of methyl Oleate, Oleic acid, Oleyl glycerols, and Oleyl glycerol mixtures in supercritical carbon dioxide, J. Am.

21、 Oil chem. Soc., 1991, 68(2): 87-91.附录 程序源代码一、模块源代码1、定义变量Option Explicit定义输入输出的文件名Public file_input As StringPublic file_weight As String定义输入和输出的数据个数Public input_number As IntegerPublic output_number As Integer定义每层的结点数Public nodes_of_layer0%, nodes_of_layer1%, nodes_of_layer2%Public learning_rate As

22、 Single输入数值和输出数值;原始数据Public inputdata() As SinglePublic outputdata() As Single归一化以后的数据Public netinput() As SinglePublic netoutput() As Single各层权重Public weight_10() As SinglePublic weight_21() As Single各层的输出值Public layer0_out() As SinglePublic layer1_out() As SinglePublic layer2_out() As SinglePublic

23、 error1() As SinglePublic error2() As SinglePublic program_state As StringPublic regulate_input As SinglePublic regulate_output As Single2、重新设置本子程序的目的:重置所有的数组,重新开始计算Sub refresh_net() Dim i%, j% ReDim inputdata(nodes_of_layer0) As Single ReDim outputdata(nodes_of_layer0) As Single ReDim netinput(node

24、s_of_layer0) As Single ReDim netoutput(nodes_of_layer0) As Single ReDim weight_10(nodes_of_layer1, nodes_of_layer0) As Single ReDim weight_21(nodes_of_layer2, nodes_of_layer1) As Single ReDim layer0_out(nodes_of_layer0) As Single ReDim layer1_out(nodes_of_layer1) As Single ReDim layer2_out(nodes_of_

25、layer2) As Single ReDim error1(nodes_of_layer1) As SingleReDim error2(nodes_of_layer2) As Single3、权重随机赋值 给各层的权重赋随机数 Randomize For i = 1 To nodes_of_layer1 For j = 1 To nodes_of_layer0 weight_10(i, j) = Rnd Next j Next i For i = 1 To nodes_of_layer2 For j = 1 To nodes_of_layer1 weight_21(i, j) = Rnd

26、Next j Next iEnd Sub4、计算正向输出结果该子程序用于正向计算输出结果Sub caculate() Dim i%, j% Dim ftemp As Single caculate layer0 For i = 1 To nodes_of_layer0 layer0_out(i) = netinput(i) Next i caculate layer1 For i = 1 To nodes_of_layer1 ftemp = 0 For j = 1 To nodes_of_layer0 ftemp = ftemp + weight_10(i, j) * layer0_out(j

27、) Next j ftemp = 1 / (1 + Exp(-ftemp) layer1_out(i) = ftemp Next i caculate layer2 For i = 1 To nodes_of_layer2 ftemp = 0 For j = 1 To nodes_of_layer1 ftemp = ftemp + weight_21(i, j) * layer1_out(j) Next j ftemp = 1 / (1 + Exp(-ftemp) layer2_out(i) = ftemp Next iEnd Sub5、调整网络权重本子程序的作用是:对网络的权重进行调整Sub

28、 adjust() Dim i%, j% Dim ftemp As Single 调整输出层的权重 For i = 1 To nodes_of_layer2 error2(i) = (netoutput(i) - layer2_out(i) * layer2_out(i) * (1 - layer2_out(i) Next i For i = 1 To nodes_of_layer2 For j = 1 To nodes_of_layer1 weight_21(i, j) = weight_21(i, j) + learning_rate * error2(i) * layer1_out(j)

29、 Next j Next i 调整隐含层的权重 For i = 1 To nodes_of_layer1 ftemp = 0 For j = 1 To nodes_of_layer2 ftemp = ftemp + weight_21(j, i) * error2(j) Next j error1(i) = ftemp * layer1_out(i) * (1 - layer1_out(i) Next i For i = 1 To nodes_of_layer1 For j = 1 To nodes_of_layer0 weight_10(i, j) = weight_10(i, j) + l

30、earning_rate * layer0_out(j) * error1(i) Next j Next iEnd Sub6、数据的保存Sub save_net() Dim i%, j% On Error GoTo 1000 Open file_weight For Output As #1 Print #1, input_number Print #1, output_number Print #1, nodes_of_layer0, nodes_of_layer1, nodes_of_layer2 Print #1, learning_rate Print #1, regulate_inp

31、ut Print #1, regulate_output For i = 1 To nodes_of_layer1 For j = 1 To nodes_of_layer0 Print #1, weight_10(i, j); Next j, i Print #1, For i = 1 To nodes_of_layer2 For j = 1 To nodes_of_layer1 Print #1, weight_21(i, j); Next j, i Close #1 Exit Sub1000: Close #1 MsgBox 无法写入配置文件,无法记录权值, vbOKOnly, Err.C

32、learEnd Sub7、自定义函数Function load_net() As String Dim i%, j% On Error GoTo 1000 Open file_weight For Input As #1 Input #1, input_number Input #1, output_number Input #1, nodes_of_layer0, nodes_of_layer1, nodes_of_layer2 Input #1, learning_rate Input #1, regulate_input Input #1, regulate_output 重新定义数组

33、refresh_net For i = 1 To nodes_of_layer1 For j = 1 To nodes_of_layer0 Input #1, weight_10(i, j) Next j, i For i = 1 To nodes_of_layer2 For j = 1 To nodes_of_layer1 Input #1, weight_21(i, j) Next j, i Close #1 load_net = OK Exit Function 1000: Close #1 MsgBox 无法打开配置文件,重新开始学习, vbOKOnly, Err.ClearEnd F

34、unction二、网络初始化及参数设置源代码Option ExplicitPrivate Sub Form_Load() Text1.Text = nodes_of_layer1 Text2.Text = Format(learning_rate, 0.00) UpDown1.Value = Text1.Text Text3.Text = input_number Text4.Text = output_number UpDown2.Value = Text3.Text UpDown3.Value = Text4.Text Text5.Text = Form1.learn_times Text

35、6.Text = Form1.max_error * 100 If Form1.method = 0 Then Option1.Value = True Else Option2.Value = TrueEnd SubPrivate Sub Image2_Click()Unload MeEnd SubPrivate Sub Image5_Click() If Text2.Text 1 Then MsgBox 学习效率的数值只能在0到1 之间，请改正, vbOKOnly + vbExclamation, 提示栏！ Text2.SetFocus End If Form1.learn_times =

36、 Text5.Text Form1.max_error = Text6.Text / 100 If Option1.Value = True Then Form1.method = 0 Else Form1.method = 1 learning_rate = Text2.Text If Text1.Text nodes_of_layer1 Or Text3.Text input_number Or _ Text4.Text output_number Then Or learning_rate Text2.Text Dim response% response = MsgBox(网络参数已被

37、改变,是否要重置网络?选定后按返回键返回。, vbYesNo + vbExclamation, 提示栏！) If response = vbYes Then Form1.total_times = 0 input_number = Text3.Text output_number = Text4.Text nodes_of_layer1 = Text1.Text 考虑到阀值的情况 If Check1.Value = 0 Then nodes_of_layer0 = Text3.Text nodes_of_layer2 = Text4.Text End If refresh_net If For

38、m1.data_ok Then save_net End If End IfEnd SubPrivate Sub Image6_Click() If Form1.data_ok = True Then refresh_net save_net End If End SubPrivate Sub Text5_Change() Option1.Value = TrueEnd SubPrivate Sub Text6_Change() Option2.Value = TrueEnd SubPrivate Sub UpDown1_Change() Text1.Text = UpDown1.ValueE

39、nd SubPrivate Sub UpDown2_Change() Text3.Text = UpDown2.ValueEnd SubPrivate Sub UpDown3_Change() Text4.Text = UpDown3.ValueEnd Sub三、学习训练窗口源代码Option ExplicitDim array1() As SingleDim array2() As SingleDim groups%Public total_times As LongPublic data_ok As BooleanPublic learn_times As LongPublic max_e

40、rror As SinglePublic method%Private Sub Form_Load() nodes_of_layer0 = 26 nodes_of_layer1 = 8 nodes_of_layer2 = 1 input_number = 1 output_number = 1 max_error = 0.01 learn_times = 10000 method = 0learning_rate = 0.9End SubFunction load_data(file_d As String, array1() As Single, array2() As Single) As

41、 Integer Dim i%, j% On Error GoTo 1000 Open file_d For Input As #1 Dim groups%, total As Long ReDim array1(1000, input_number) ReDim array2(1000, output_number) groups = 1 total = 0 regulate_input = 0 regulate_output = 0 Do While Not EOF(1) For i = 1 To input_number Input #1, array1(groups, i) If re

42、gulate_input Abs(array1(groups, i) Then regulate_input = Abs(array1(groups, i) total = total + 1 Next i For j = 1 To output_number Input #1, array2(groups, j) If regulate_output Abs(array2(groups, j) Then regulate_output = Abs(array2(groups, j) total = total + 1 Next j groups = groups + 1 Loop1000:

43、If EOF(1) Then regulate_input = regulate_input * 2 regulate_output = regulate_output * 2 groups = total / (input_number + output_number) load_data = groups Close #1 Exit Function Else Close #1 MsgBox 无法打开数据文件, vbOKOnly, Err.Clear End IfEnd FunctionPrivate Sub Image1_Click() If data_ok Then save_netE

44、nd SubPrivate Sub Image2_Click() Dim i%, j% Dim x% If Not data_ok Then MsgBox 请先打开数据文件!, vbOKOnly, Exit Sub End If Dim ftemp As Single Select Case method Case 0: For x = 1 To learn_times ftemp = 0 For i = 1 To groups For j = 1 To input_number netinput(j) = array1(i, j) / regulate_input Next j For j

45、= 1 To output_number netoutput(j) = array2(i, j) / regulate_output Next j caculate For j = 1 To output_number ftemp = ftemp + Abs(layer2_out(j) - netoutput(j) / netoutput(j) Next j adjust Next i DoEvents Image5.Visible = True total_times = total_times + 1 Text2.Text = total_times ftemp = ftemp / (ou

46、tput_number * groups) Text1.Text = Format(ftemp * 100, 0.000000000) & % Next x Case 1: Do ftemp = 0 For i = 1 To groups For j = 1 To input_number netinput(j) = array1(i, j) / regulate_input Next j For j = 1 To output_number netoutput(j) = array2(i, j) / regulate_output Next j caculate For j = 1 To output_number ftemp = ftemp + Abs(layer2_out(j) - netoutput(j) / netoutput(j) Next j adjust Next i ftemp = ftemp / (output_number * groups) Text2.Text = tota