MATLAB在相平衡计算中的应用

曾小红

(江西省化学工业设计院，江西南昌 330003)

近年来，国内一些研究人员开始尝试使用MATLAB最优化工具箱(optimization toolbox)来研究VLE。苏铁军［1，2］运用MATLAB最优化工具箱中的非线性最小二乘法程序，对二元和三元汽液平衡进行Wilson模型参数估算。计算结果表明这种方法大大减轻了复杂算法设计的工作量。刘俏［3，4］用MATLAB计算化工过程中的汽液平衡组成及温度，结果表明MATLAB具有编程容易、计算方便和效率高的特点。唐正姣等［5］运用MATLAB最优化工具箱函数对二元系等温等压汽液平衡数据进行UNIQUAC模型参数估算。

本文选用101.325kPa下的甲基二氯硅烷(1)+甲苯(2)［6］二元体系为研究对象，用MATLAB最优化工具箱中的非线性最小二乘法来进行Wilson模型参数回归和泡点计算。

1.非线性最小二乘法及其函数

非线性最小二乘法的MATLAB求解函数为lsqnonlin()。非线性最小二乘法在数据拟合和参数估计等方面有广泛用途。

数学模型:

lb≤x≤ub式中，C为常数。

调用格式:

输入函数:

fun 函数返回值fun(x)

x0 初值

lb，ub 设置变量的下边界和上边界，lb≤x≤ub，若无边界，则将lb和ub赋给空矩阵:lb=［］，ub=［］。如果x(i)下无界，则设置lb=－inf;如果x(i)上无界，则设置ub=inf;

输出参数:

x 返回解向量;

resnorm 返回x处残差的平方范数值:sum(fun(x).^2);

residual 返回x处的残差值fun(x);

exitflag 描述函数计算的退出条件:若为正值，表示目标函数收敛于解x处;若为负值，表示目标函数不收敛;若为零值，表示已经达到函数评价或迭代的最大次数。

output 返回优化信息:

output.iterations 表示迭代次数;

output.algorithm 表示所采用的算法;

output.funcCount 表示函数评价次数。

Lambda 返回包含x处拉格朗日乘子的结构参数;

jacobian 返回解x处的fun函数的雅可比矩阵。

lsqnonlin默认时选择大型优化算法。lsqnonlin通过options.LargeScale设置为‘off’来作中型优化算法。其采用一维搜索法。

2.程序说明

本文采用 MATLAB7.0 编程，MATLAB［7］具有一系列专题应用软件包，称之为工具箱(Toolbox)。工具箱是极为广泛的集合了MATLAB函数(M－files)，而这些函数大大地扩展了MATLAB环境，使之能够解决某些专门的问题。工具箱所涉及的领域，包括信号处理、控制系统设计、鲁棒控制、动态系统仿真、系统辨识、神经网络、图像处理、统计学、μ分析、模型逻辑、优化、样条函数、符号运算、图形设计、符号数学处理、偏微分方程处理、高阶谱分析等。

本文利用MATLAB语言提供的最优化工具箱可以直接调用lsqnonlin函数模块，因而使得程序的编写和调试非常简便。本文编程中用到的命令如表1－1所示:

表1－1 程序中所用的基本函数表

3.Wilson模型参数回归

低压下非理想溶液(汽相看做是理想气体)的汽液平衡关系式:

纯组分i的饱和蒸气压Psi(单位为mmHg)由Antoine方程计算:

式中 A、B、C ——Antoine参数。

Wilson模型活度系数方程:

纯液体摩尔体积 VL1，VL2由修正的 Rackeet方程［8］计算:

拟合程序的目标函数可以用多种方式来定义。拟采用以下函数φ对液相活度系数的Wilson模型参数进行参数回归:

式中Pcal和Pexp分别表示混合物中某组分在平衡温度下蒸汽压的计算值和实验值。

4.泡点计算

泡点计算主要求取液相混合物在一定压力下的沸点即泡点，或一定温度下的蒸汽压，以及平衡汽相组成。输入的独立变量为x以及T或P的任一个，要求输出另一个以及y。这些问题都是单变量非线性方程的求解问题，本文用MATLAB的主程序fzero()求解，并以温度实验值作为其计算的初值。

等压泡点计算程序见图4.1:

图4.1 泡点温度与汽相组成的计算框图

5.结果与讨论

5.1 纯物质的Antoine常数及临界性质

从文献中查得的Antoine参数及临界参数示于表5－1。

表5－1 试剂Antoine方程参数及临界参数:lg(P/KPa)=A－B/［(T/K)+C］

5.2 MATLAB最小二乘法Wilson模型参数回归结果

表5－2 模型参数

5.3 模型误差比较

表5－3 各种模型计算二元体系VLE的误差

5.4 二元体系的VLE数据

最小二乘法计算二元系的汽液平衡数据关联的结果示于附表5－1，从表5－3中可看出，最小二乘法的平均绝对偏差小于文献计算的结果，

附表5－1 101.3kPa下甲基二氯硅烷－甲苯二元系的实验值和Wilson模型计算值

No. x1exp Texp Tmod dT y1exp y1mod dy 24 0.6188 329.91 329.72 0.19 0.9443 0.9461－0.0018 25 0.6485 328.40 328.18 0.22 0.9523 0.9541－0.0018 26 0.6812 326.89 326.54 0.35 0.9598 0.9619－0.0021 27 0.7146 325.29 324.94 0.35 0.9671 0.9689－0.0018 28 0.7467 324.02 323.46 0.56 0.9724 0.9747－0.0023 29 0.7820 322.48 321.91 0.58 0.9784 0.9803－0.0019 30 0.8137 321.27 320.57 0.70 0.9827 0.9846－0.0019 31 0.8429 320.05 319.39 0.66 0.9867 0.9881－0.0014 32 0.8733 318.80 318.23 0.57 0.9904 0.9912－0.0008 33 0.9057 317.72 317.05 0.67 0.9933 0.9941－0.0008 34 0.9385 316.67 315.92 0.75 0.9959 0.9965－0.0006 35 0.9651 315.85 315.07 0.78 0.9978 0.9982－0.0004 36 0.9830 315.31 314.53 0.78 0.9990 0.9992－0.0002 Max 1.27 0.0294 MAE 0.55 0.0056

［1］李克娥，苏铁军.MATLAB在处理化工热力学数据中的应用［J］.重庆科技学院学报(自然科学版)，2005，7(3):100 －102.

［2］苏铁军.利用MATLAB计算三元物系的汽液平衡数据［J］.江汉大学学报(自然科学版)，2005，33(4):32－34.

［3］刘俏，范圣第.基于 MATLAB 的化工计算［J］.计算机应用与软件，2005，22(12):130 －132.

［4］刘俏.MATLAB辅助化工原理教学的时间［J］.高等理科教育，2006，5:80 －83.

［5］唐正姣，王存文，柯文语，等.用MATLAB估算二元系的 UNIQUAC模型参数［J］。2007，45(3):8－10.

［6］Yu SX，Yu ZB，Qiu ZM，et al.Isobaric Vapor—liquid Equilibrium for Methyldichlorosilane－methyl Vinyldi-chlorosilane－toluene and Constituent Binary Systems［J］.Chinese Journal of Chemical and Engineering，2003，11(2):213 －216.

［7］飞思科技产品研发中心编著.MATLAB6.5辅助神经网络分析与设计［M］.北京:电子工业出版社，2003:1－12.

［8］Harold G Rackett.Equation of state for saturated liquids［J］.Journal of Chemical and Engineering Data，1970，37:514 －517.

［9］谢心亮.多元卤代硅烷体系的汽液相平衡研究［D］.南昌大学，2005.