D0I:10.13374/i.issn1001-053x.1992.05.002 第14卷第5期 北京科技大学学报 Vol.14No.5 1992年9月 Journal of University of Science and Technology Beijing Sept.1992 最小自由能原理的SUMT方法° 郭汉杰· 赵玉祥· 精要:用改进的SUMT方法编写的FORTRON程序,对等式约束的化学平衡的最小自由能问题 求解,精度优于同类程序,并可用于普通微机。改进的主要部分是强迫等式约束达到一定的 精度后,主程序进人收敛判断, 关键词:SUMT,最小自由能,化学平衡 Method of Minimizing the Global Free Energy with SUMT Program Guo Hanjie· Zhao Yuriang· ABSTRACT:By improving on the FORTRON program of SUMT,the algorithm of minimizing the global free energy of the system (G),which is the nonlinear programming with equality constraints, have solved.The precision of the algorithm is higher than the others.The program suit to PC comput- ers.The main parts improved on the program is that equality constraints must be forced to a specified precision,before the main program converges. KEY WORDS:SUMT,minimizing the global G,chemical equilibrium 多相多组分体系的化学平衡的计算已有不少文章发表”,而直接利用非线性规划方法求 算平衡组分可以避免VC$2”方法中对初值的要求苛刻、计算时间长等缺点。国外已有一些用非 线性规划方法求算化学平衡的程序),但由于技术保密,其原程序不易得到。所以开发一种适 用于微机的求算多相多组分体系的化学平衡的通用程序是必需的。 多相多组分体系的化学平衡的最小自由能原理的数学模型是: 工逃) min.G=22xn(g+RTln k=1j=1 =】 ①1992-01一20牧稿 冶金系(Department of Metallurgy) ·502·
第 卷第 期 年 月 北 京 科 址 技 大 学 学 报 恤 最小 自由能原理 的 方法 ① 郭汉 杰 ‘ 赵玉 祥 ’ 摘要 用改进的 方法编写 的 程 序 , 对等式 约束 的化学平衡 的最 小 自由能 问题 求解 , 精度优于同类程序 , 并可用于普通微机 。 改进的主要部分是强迫等式约 束达到 一定的 精度后 , 主程 序进人收敛判 断 。 关徽词 , 最小 自由能 , 化学平衡 召别。 ,绝 ‘ 么初。 触劣心 ’ 加讲 讲。歹 , , , 少 七 琳到治 写 , 刀旧 , 妞址协 , 多相 多组分体系的化学平衡的计算已有不少文章发表 〔一 〕, 而直接利用 非线性规划方法求 算平衡组分可以避免 “②方法 中对初值的要求苛刻 、 计算时 间长等缺点 。 国外 已有一些用非 线性规划方法求算化学平衡的程序 〔的 , 但 由于技术保密 , 其原程序不 易得到 。 所 以开发一种适 用于微机 的求算多相多组分体系的化学平衡的通用 程序是必需 的 。 多相 多组分体系的化学平衡的最小 自由能原理 的数学模型是 。 一 全云。 衅 井 仓二 夕 刃 二产 ① 一 一 收稿 , 冶金系 众,功 堆 · · DOI:10.13374/j.issn1001-053x.1992.05.002
s.t.h=)-b.=0 (1) k=13=1 i=1,2…,C x≥0〔(j=1,2…,),(k=1,2,…,m)〕 这是一个非线性规划问题(1)。目前解此类问题比较好的软件有POPI、NLP、GGS、GRG、 SUMT及可变容差法?,但前5个软件需到美国有关公司去买,能得到的软件有可变容法和国 内公开的SUMT)等。 1用现有的软件对问题(1)的计算结果 对如下反应平衡: H20+C0=H2+C02 (2) 问题(1)中参数值为:m=1,=4,C=3,即单相、·4组元、3元素的体系,900K时: 唱2=-128.80kJ 8o=-299.84kJ 20=-423.89kJ 80,=-602.25kJ) 由化学平衡计算得,H20、C0、H2、C02平衡组分分别为:(初始H2O,C0的量为1mol) H,0=co=0.4087mol 2=xco2=0.5913mol 用文献〔6)公布的可变容差程序计算结果为: ,2193451E-04 6318092E-05 .4482269E-04 THESE ARE FINAL ANSWERS ¥(1)=.592365E+00 x(2)=.407646E+00 x(3)=.407614E+00 x(4)=.592392E+00 Stop-Program terminated. 其中x(1),x(2),x(3),x(4)分别为H0、C0、H2、CO2的平衡成分。上面3个数 字,2193451e等为计算精度,即3个等式约束的精度。可以看出可变容差法计算值与化学 平衡法计算值非常一致。 用文献〔T)公布的SUMT程序的计算结果如下:其中上面为输入的初值,变量上、下限、 精度、、及3个等式约束的系数矩阵所组成的数据文件。结果表中x(1)~x(4)分别表示 H2、Hz0、CO、C02的平衡摩尔数。 x(1)=.242388E+00 9(1)=.24388 x(2)=.134187E+00 g(2)=1.341867 x(3)=.234510E+00 9(3)=.234510 x(4)=.194829E+01 g(4)=1.948288 f(1)=一250.955500 h(1)=.116851E+01 ·503·
· 二 刃 刃 凡声 , 一 阮 正 一 , … , 、 , 〔 , … , , 一 , , … , , 〕 这是一个非线性规划 问题 。 目前解此类 问题 比较好的软件有 、 、 、 、 及可变容差法 ‘幻 , 但前 个软件需到美国 有关公 司去买 , 能得到 的软件有可变容法和 国 内公开 的 〔,等 。 用现有的软件对问题 的计算结果 对如下反应平衡 问题 中参数值为 二 一 , 一 , , 即单相 、 组元 、 元素的体系 , 时 嵘 一 · 碍 。 一 一 · 嵘 二 一 · 。 一 〔“〕 由化学平衡计算得 , 、 、 、 平衡组分分别为 初始 , 的量为 ‘ 一 鲡 一 · ‘, 」 用文献 〔印 公布的可变容差程序计算结果为 一 一 犷 一 十 , 一 一 其中 , , , 分别为 、 、 、 的平衡成分 。 上面 个数 字 一 ‘等为计算精度 , 即 个等式约 束的精度 。 可以看 出可变容差法计算值与化学 平衡法计算值非常一致 。 用 文献 〔〕公布的 程序的计算结果如下 其 中上面 为输入的初值 , 变量上 、 下限 、 精度 、 衅 、 及 个等式约束的系数矩 阵所组成的数据文件 。 结果表 中 分别表示 、 、 、 的平衡摩尔数 。 十 夕 二 十 艺 夕 劣 夕 一 十 · ·
h(2)=.118280E+01 h(3)=.347295E+01 看来国内的SUMT程序对处理具有等式约束的非线性规划问题还很不成熟。可变容差法在变 量数增多时计算时间迅速增长,如本文第3部分提出的问题中,可变容差法在P℃/XT微机上 运算10h也未收敛。显然这种方法是不可取的。 2改进的SUMT方法 SUMT(Sequential Unconstrained Minimization Technique)是序列无约束极小化方法的简写, 对问题(3): min.f() x∈B s.t.9:(2)≥0 i=1,2,9 (3) h,()=0 j=1,2,…,p(p<n) 构造的惩罚函数的具体形式为: Pw)-+名高+元三 (4) 式中为罚因子,是一个递减的无穷正数数列。若r=,极小化罚函数P(,),可求得 相应的极值点x(r)。对一序列罚因子{},当k→∞时,→0,在适当条件下,其相应 的极小点x()使: =蓝高=0 (5) ={高r=0 P (6) 即有: mP(i(r)=f(在) (7) 求罚函数P(i,r)的无约束极小化采用Pow©l的搜素法cm,,这是一种不用导数的无约束极 小化方法,它比导数法有明显的优点。因为导数法在维数稍高时,为梯度法和二阶导数法提 供解析函数是困难的或者是不可能的。另外Powell法不要求目标函数的正则性、连续性和导 数的存在性。对问题(3),SUMT调用Powell的基本方法是,从初始点x)沿着由过程生成的 一组共轭方向作逐次一维搜索。在逐次产生的搜索方向组忌,(=1,2,…,)中,为了避免 个搜索方向的相关性,每次求得共轭方向京+1后,不等式: f3<f (8) 与 (f一2f2+f3)(f-f2-4)2<0.54(f1-fs)2 (9) 同时成立,则用3+1方向代替3。方向,否则仍用原方向组进行迭代求解。式中1为初值点x0 的函数值;∫2为沿3.方向搜索点云.的函数值;f为沿共轭方向京+1上反射点+1(云+1=2云 一xo)的函数值:4为及(i=1,2,…,)方向搜索中最大的函数值差,即: ·504·
十 看来国 内的 程序对处理具有等式约 束的非线性规划 问题还很 不成熟 。 可变容差法在变 量数增 多时计算时间迅速增长 , 如本文第 部分提出的 向题 中 , 可变容差法在 微机上 运算 也未收敛 。 显然这种方法是不可取的 。 改进的 方法 对问题 们口们 。 让 是序列无约 束极小化方法的简写 , 劣 ,‘王 妥任 矛 , , … , 凡 构造 的惩罚 函数的具体形式为 “ , , … , 尹 尹 论 今 , 劣 ,, , 刃 ,万,戈… 十 二二二 ‘、 少 了 刀 〔无, 〕 式 中 , 为罚因子 , 相 应 的极值点 妥 是一个递减的无穷正数数列 。 若 , 一 , , 极 小化罚 函数 ‘ , , , 可求得 , 。 对一序列罚因子 伽 , 当 时 , 产 。 , 在适 当条件下 , 其相应 的极 小点 奋 沪 , 使 、、产 了飞 勺内︻ 尸、、 口 枷 , 互,刃 盈 , ‘二 架‘两 刃 〔凡 〕 一 即有 忽妥,“,一 了王 求罚 函 数 ‘ , ,‘ 的无约束极小化采用 的搜索法⑦ , 这是一种不用 导数的无约束极 小化方法 , 它 比导数法有明显的优点 。 因为导数法在维数稍高时 , 为梯度法和 二 阶导 数法提 供解析 函数是困难的或者是不可能的 。 另外 法不要求 目标 函数的 正则性 、 连续性和导 数的存在性 。 对问题 , 叮 调用 的基本方法是 , 从初始 点 扮 , 沿着 由过程生成的 一组共扼方 向作逐次一维搜索 。 在逐次产生的搜索方 向组 瓦 一 , , … , 二 中 , 为了避免 。 个搜索方 向的相关性 , 每次求得共扼方向 反 后 , 不等式 九 与 一 、 一 一 , 二 、 一 了 同时成立 , 则用 反 方向代替 若 。 方 向 , 否则 仍用原方 向组进行迭代求解 。 式 中 了 , 为初值点 乳 的函 数值 介 为沿 瓦方向搜索点 王 的 函数值 介 为沿共扼方 向 反 十 上反射点 牙、 , 三 二 云 一 乳 的 函数值 」 为 及 ‘ , , 一 , 方 向搜索中最大的函数值差 , 即 · ·
4='{}=恶{P()-P(-1)} (9a) 京。为与4相对应的方向。 无约束极小化的收敛准则为: fz-x1)月≤e (10) 若还不是约束最优点,则按惩罚函数法用新的,-D构造罚函数P(在,+D),并求其无约 束极值点,这样便产生了一系列的点列: x(r),z(2),…,(rD,z(r) 当罚函数值满足 P((r)-P((r-) P(x(r-)) ≤z (11) 时,SUMT已收敛,认为:(r)点即为问题(3)的约束最优解。图1列出了SUMT调Powel 的程序框图,用该程序计算本文第1节所提出的问题显然是不行的,对仅有不等式约束的非 线性最优化本程序有一定的精度,而对于具有等式约束的非线性规划问题有必要进行改进。 为了使图1适应于等式的约束问题,在图1()~(b)段插入图2段,该段功能是强迫 等式约束达到如下精度 H.= |H()<e1 (12) 的前提下,主程序才做收敛判断式(11)。 该程序段的加入收到了满意的效果。重新计算问题(2),结果如下: x(1)=,591294E+00 g(2)=.408706 x(2)=.408706E+00 g(3)=.408704 x(3)=.408704E+00 g(4)=.591296 x(4)=.591296E+01 h(1)÷.000000E+00 f(1)=-100.236600 h(2)=-,834465E-06 g(1)=.591294 a(3)=.953674E-06 与平衡常数法所得结果完全一致。 3改进的SUMT与著名软件的比较 Himmelblau)曾调用了美国几个著名的非线性规划软件:NLP、可变容差法、GGS、GRG、 SUMT求算如下问题(等温等压下化学平衡): 独立变量数:10 约束:线性等式的约束3个 目标函数: 10 10 fx〕=zx(C十ln(x/zx)) (3) C1=-6:089 C2=-17.164 ·505·
一 血 妄 ‘一 尸妄卜 吸卜屯 瓦 为与 』 相对应的方向 。 无约束极小化的收敛准则为 忏 一 会卜 月簇 街 若 补 ,还不是约 束最优点 , 则按惩罚 函数法甩新的 ,卜‘,构造罚 函数 孰 ,歼 ‘ , 束极值点 , 这样便产生 了一系列 的点列 妥,“, , 妥, ‘,, , … , 主,‘,一 ‘,,奋,“, 当罚函数值满足 ‘ 并求其无约 尸妥,‘,,一 尸奋‘,一 ‘, 公 护卜‘ 、 。 时 , 已收敛 , 认为 妥秘 点 即为问题 的约束最优解 。 图 列出了 调 的程序框图 , 用该程序计算本文第 节所提 出的间题显然是不行的 , 对仅有不等式约束 的非 线性最优化本程序有一定的精度 , 而对于具有等式约束的非线性规划 问题有必要进行改进 。 为 了使 图 适应于等式的约束间题 , 在 图 段插入图 段 , 该段功能是强迫 等式约束达到如下精度 拒万石 的前提下 , 主程序才做收敛判断式 。 该程序段的加人收到 了满意 的效果 。 劣 劣 劣 一 重新计算问题 , 结果如下 , 夕 夕 十 一 一 一 与平衡常数法所得结果完全一致 。 改进的 与著名软件的比较 珍加 扭 曾调用了美国几个著名的非线性规划软件 、 可变容差法 、 、 求算如卞间题 等温等压下化学平衡 独立变量数 约束 线性等式 的约束 个 目标函数 , 、 〔〕二 刀 ‘‘ 垃 二‘刀 , ‘二 二 二 一 一 一 · 辉产 训 洲
Star nKr.Kg.Km.b(k).bu(k)x(k).cnez.r.c.f x(使)=x(k).m(k)=x(k) Cy=Cy+1 456 Call "funct",compute f9()( 〔k)<1015 Call "Powe 11",computex(n(kg Call "RANDOM",determine x(kk=1…刀 9(k10-=上 Zm(x( Print the symbol"" ..."of abuormal end 122 012 Print the symbol "[ .."of normal fnd end b rar+C xm(k)=x(k) 125 图1SUMT调Powel原始程序框图 图2在原SUMT增加的程序段 Fig.1 The original Schematic diagram of SUMT calling"Powell" Fig.2 The adding program on the original program C3=-34.054 C4=-5.914 C5=-24.721 C6=-14.986 C,=-24.100 Cg=-10.708 Cg=-26.662 C10=-22.179 约束: a1〔)=x1十2x2+2xg+x6+10-2=0 2[=x4十2x5十6+x1-1=0 ·506·
﹁‘,备 , 了仍 知雨 ‘二二 目 厂,﹄ 图 调 原始程序框 图 图 在原 增加 的程序段 堪如 乎 】 即 开 招 琢山 胡 一 一 二 一 口。 一 二 一 。 一 , 二 一 一 约束 无〔易〕 二 。 。 一 一 。 、〔妄〕 二。 。 二 一 二 ·
3Cr)=x3+x7+xg+2xg十x10一1=0 x≥0 元=1,2,,10 不可行起始点: x0)=0.1 元=1,2,w…,10 用本文提供的改进的SUMT计算结果见表1。 可以看出,改进的SUMT程序精度已优于国外同类程序,最优解已完全收敛于等式约束 的线上。在PC/XT微机上运算时间只需5min即可收敛。 表1对问题(3)美国几个著名程序的计算结果 Table 1 The results of problem (3)computed by.some famous programs in USA NLP 可变容差法 GGS GRG SUMT f(] -47.751 -47.736 -47.656.,-47.761 -47.761. 1 0.0350 0.0128 0 0.0406 0.0407 0.1142 0.1433 0.1695 0.1477 0.1477 g 0.8306 0.8078 0.7536 0.7832 0.7832 4 0.0012 0.0062 0. 4 0.0014 0.0014 0.4887 0.4790 0.5000 0.4853 0.4853 、 6 0.0005 0.0033 0 0.0007 0.0007 7 0.0209 0.0324 0 0.0274 0.02747 0.0157 0.0281 0 0.0180 0.0180. 珀 0.0289 0.0250 0.0464 0.0375 0.0373 x10 0.0751 0.0817 0.1536 0.0969 0.0969 h(r) 3.E-12. 3.E-05 0 1.E-06 -8.E-08 h2〔x 3.E-12 2.E-05 0 1.E-06 -1.E-07 (x) 2.E-11 9.E-05 0 1.E-06 -1.E-07 另外,用本程序对20个变量和30个变量的问题分别进行试算,占用机时分别为100min 和5h。 x(1)=.345405E-01 9(2)=.145416 x(2)=.145416E+00 g(3)=.787161 x(3)=.787161E+00 g(4)=.004122 x(4)=.412205E-02 9(5)=,487011 x(5)=,487011E+00 g(6)=.001943 x(6)=.194348E-02 g(7)=.019913 x(7)=.199132E一01 g(8)=.014951 x(8)=.149511E-01 9(9)=.039806 x(9)=.398063E-01 g(10)=.098362 x(10)=.983621E-01 h(1)=.000000E+00 f(1)=-47.756640 (2)=.000000E+00 g(1)=.034541 h(3)=.000000E+00 ·507
。〔〕 公 。 。 。 一 一 坛 , , … , 不可行起始点 对 ‘一 , , 。 , 用本文提供的改进的 计算结果见表 。 可 以看出 , 改进 的 程序精度 已优于 国外同类程序 , 最优解已完全收敛于等式约束 的线上 。 在 微机上运算时 间只需 两 即可收敛 。 ‘ ’ 表 对问题 美国几个著名程序的计算结果 , ,功 熟功碌笋 卜 积 协 可变容差法 一 ‘ 了〕 二 公 忿 雷 劣 工 忿, 工 工 〔 劣〕 〔〕 几〔 子 〕 一 一 。 。 。 ‘ 。 ‘ , 一 一 。 一 。 一 。 一 一 一 一 。 。 盯 一 魂 · 。 ‘ 石日 一 一 。 一 一 ‘ 、 了 、 一 吐 。 。 公 勺 一 一 一 一 一 , 一 另外 , 用本程序对 个变量和 个变量的问题分别进行试算 , 占用 机时分别为 和 。 劣 一 二 十 夕 二 一 一 夕 一 公 二 一 夕 一 夕 二 一 一 一 一 夕 二 · ·
4结论 改进的SUMT方法,对用最小自由能原理求解化学平衡是行之有效的,计算精度已优于 国外同类程序;在微机上运算结果表明,30个变量以下的化学平衡问题完全可以在微机上运 行,10个变量在Pc/XT机上占机时5min,20个变量约100min,30个变量5h。 附:全文符号: G:体系总自由能;m:相数;:组元数,C:元素数;x:无相中j组元摩尔数:g:k相 中j组元的标准摩尔自由能;R:气体常数;T:绝对温度(K);五(x):元素的等式约束向 量;E:k相中i元素在j组元中的系数;b,:i元素的总量;9:():不等式约束的表达式向 量;T:惩罚因子;P(x,r):惩罚函数;f(2):目标函数;e1,2,e1:分别为收敛准则规定 的精度,一般为104~105。 致谢:忠谢何晓峰博士对本工作所提供的帮助! 参考文献 .1 Smith W R.Theoretical Chemistr Advance and Perspectives.New York:Academic Press, 1980,5 2 Villars D S.J of Physical Chemistry,1959,63:521 3许志宏,王乐珊.无机热化学数据库.北京:科学出版社,1987 4 Gautam C W.AICHE J,1979,25:991 5 White C W.AICHE J,1981,27:466 6 Himmelblau D M.Applied Nonlinear Programming.New York:McGraw Hill Book Compa- ny,1972 7万耀青,梁庚荣.最优化计算方法常用程序汇编.北京:工人出版社,1983 8 Burin I,Knacke O.Thermochemical Properties of Inorganic Substances.USA:Springer Press 73),Supplement ('77 )1977 ·508·
结论 改进的 方法 , 对用最小 自由能原理求解化学平衡是行之有效的 , 计算精度已优于 国外 同类程序 在微机上运算结果表 明 , 个变量 以下的化学平衡问题完全可 以在微机上运 行 , 个变量在 机上 占机时 , 个变量约 , 个变量 。 附 全文符号 体系总 自由能 相数 , 组元数 元素数 孙 相 中 了组元摩尔数 嵘 相 中 组元的标准摩尔 自由能 气体常数 绝对温度 , 勃 ‘元素的等式约束 向 量 , 相中 元素在 ,组元 中的系数 ‘ ‘元素的总量 , 幻 不等式约束的表达式 向 量 , 惩罚 因子 尸 ‘ , , 惩罚 函数 幻 目标函数 。 , 。 , 。 分别为收敛准则规定 的精度 , 一般为 一。 一 ‘一 一 。 致谢 感谢何晓峰博士 对本工作所提供的帮助 参 考 文 献 】 , , 辽 , , 许志宏 , 王乐珊 无机热化学数据库 北京 科学 出版社 , 饭 , , 锐 , , 如恤 、 凡 致阅 , 万耀青 , 梁庚荣 最优化计算方法常用程序汇编 北京 工人 出版社 , , 而份川 比 ‘ , ’ ,