SAT-based Bounded and Unbounded Model Checking

第20卷第6期2020年12月交通运输系统工程与信息Journal of Transportation Systems Engineering and Information TechnologyV ol.20No.6December 2020文章编号：1009-6744(2020)06-0170-08中图分类号：U268.6文献标志码：ADOI:10.16097/ki.1009-6744.2020.06.022城市轨道交通列车追踪间隔与牵引能耗优化高豪，郭进*，张亚东(西南交通大学信息科学与技术学院，成都611756)摘要：针对城市轨道交通高峰小时列车密集追踪运行的特点，将降低列车牵引能耗和提升线路通过能力同时作为优化目标，研究列车运行操纵优化问题.给出移动闭塞条件下列车牵引能耗和最小追踪间隔的计算方式，考虑列车安全、正点运行约束，构建双目标优化模型.结合ε-约束法，提出一种基于动态规划的搜索算法求解模型.以亦庄线为优化算例，求解得到一组列车最优操纵Pareto 解，体现两优化目标之间的均衡关系：列车进站过程采用两次制动操纵策略可有效压缩最小追踪间隔，为弥补两次制动过程额外消耗的运行时间，列车需付出更多的牵引能耗提升进站以前的运行速度以满足正点运行约束.关键词：铁路运输；列车最优操纵；动态规划；列车；最小追踪间隔；节能Optimization of Train Headway and Traction EnergyConsumption in UrbanRail TransitGAO Hao,GUO Jin,ZHANG Ya-dong(School of Information Science and Technology,Southwest Jiaotong University,Chengdu 611756,China)Abstract:Focusing on the characteristics of dense train tracking operation during the peak hour in urban rail transit,this paper takes both reducing traction energy consumption and improving line capacity as objectives in the train driving strategy optimization.The calculation method of traction energy consumption and minimum headway is given under the moving block system.Considering the constraint of safe running and punctual arriving,the two-objective optimization model is constructed.A dynamic programming based approach combined with ε-constraint method is proposed to solve the model.A case study of Yizhuang urban rail line was executed and a set of Pareto solutions were achieved which reflecting the trade-off between those two objectives.The results showed that the multi-step braking strategy during the phase of entering a station can compress the minimum headway effectively.To reduce the increased running time because of the multi-step braking strategy,trains need to improve their velocity before the phase of entering a station to satisfy the constraint of punctual arriving,which consumes more traction energy.Keywords:railway transportation;optimum driving strategy;dynamic programming;train;minimum headway;energy saving0引言随着城市轨道交通大量新建线路投入网络化运营，我国轨道交通客运量增长明显，北京、上海两市日均客运量均已超过1000万人.为满足日益增长的客运需求，以基于通信的列车运行控制(Communication Based Train Control,CBTC)为代表的移动闭塞系统取代了传统固定闭塞系统，并得到广泛应用.利用高精度列车定位、双向大容量车地无线通信等新兴技术，CBTC 系统可以进一步提升高峰小时线路通过能力.据资料统计，2019年收稿日期：2020-07-15修回日期：2020-09-16录用日期：2020-09-17基金项目：国家自然科学基金青年科学基金/Young Scientists Fund of the National Natural Science Foundation of China(61703349).作者简介：高豪(1989-)，男，江苏苏州人，博士生.*通信作者：*******************.cn第20卷第6期城市轨道交通列车追踪间隔与牵引能耗优化全国城市轨道交通高峰小时最小发车间隔平均为290s，进入120s以内的线路共12条.同线路能力一样，列车牵引能耗是制约城市轨道交通可持续发展的关键因素.面向高峰小时多列车密集追踪的运行场景，能耗问题尤为突出.在保证列车安全、正点运行前提下，通过调整列车操纵策略，进一步压缩行车间隔并减少牵引能耗，可有效提升城市轨道交通服务水平和经济效益.针对多列车运行操纵优化问题，国内外学者开展了广泛讨论.Wang[1]面向固定闭塞和移动闭塞系统下的列车追踪运行场景，以最小追踪间隔为约束条件，研究列车运行受扰后的多列车节能操纵优化问题.Ye[2]考虑快慢车越行场景下的列车安全追踪及正点运行约束，通过构建多阶段优化控制模型求解多列车节能操纵最优策略.进一步地，柏赟[3]考虑追踪间隔要求和再生制动能利用，研究以列车净能耗最少为目标的快慢车线路列车协同操纵优化问题.Wang[4]在确保列车运行满足最小追踪间隔的前提下，通过调整各列车的站间运行时间来最大化多列车追踪优化的节能效果.以上文献将最小追踪间隔作为约束条件，研究不同场景下多列车节能驾驶优化问题.优化后的列车运行间隔发生改变，但线路能力并未得到本质提升.Takeuchi[5]将最小追踪间隔作为衡量线路能力的性能指标，分析列车运行速度对最小追踪间隔的影响.陈荣武[6]通过调整车站限速区域及限速值压缩了近6.3%的列车最小追踪间隔.Nakamura[7]指出列车采用多级制动的进站操纵策略可有效压缩最小追踪间隔，并利用遗传算法优化列车进站操纵策略.上述文献将最小追踪间隔作为优化目标，通过调整列车操纵策略提升线路通过能力，却忽略了操纵策略改变带来的牵引能耗增加问题.针对高峰小时地铁列车密集追踪运行的特点，本文同时将降低列车牵引能耗和缩短最小追踪间隔作为优化目标，通过优化列车操纵策略实现列车节能驾驶并提升线路通过能力.首先给出移动闭塞条件下列车牵引能耗和最小追踪间隔的计算方式，在此基础上考虑列车正点运行和线路限速约束，构建双目标优化模型；引入ε-约束法将模型转为单目标形式，并利用动态规划方法作进一步求解.1列车运行性能指标1.1列车牵引能耗对运行过程中的列车进行受力分析，构建基于单质点的列车运动学计算模型，即d v()xd x=u()x-r()xM⋅()1+α⋅v()x(1)r()x=r b()x+r1()x(2) r b()x=M⋅[]a+b⋅v()x+c⋅v2()x(3)r1()x=r g()x+r c()x+r t()x(4)式中：M为列车质量；x为列车位置；v()x为列车速度；α为回转系数；u()x为列车牵引制动力；r()x为列车运行阻力，包括基本阻力r b()x和和附加阻力r1()x；r b()x由戴维斯方程来表示，其系数a、b及c根据列车型号而定；r1()x包括坡道附加阻力r g()x、曲线附加阻力r c()x和隧道附加阻力r t()x.列车以操纵策略U={u()x|u min()x≤u()x≤}u max()x,x∈[0,X]在线路上追踪运行，其中，u max()x 和u min()x分别为列车保证乘客安全、舒适条件下的最大牵引力和最大制动力，X为列车站间运行距离.列车执行U后的牵引能耗E()U和运行时间T()U为E()U=∫0X max[]0,u()x⋅d x(5)T()U=∫0X1v()x⋅d x(6) 1.2最小追踪间隔线路通过能力体现为列车追踪运行过程中最小追踪间隔的倒数[5]，取决于列车在线路最受限制点处的最小安全追踪距离及通过该距离的运行速度[6].移动闭塞系统中，相邻列车的间隔距离必须始终大于最小安全追踪距离.现有CBTC系统均采用不考虑先行列车速度的“硬撞墙”模型来计算最小安全追踪距离，即S min()x a=L react()x a+L eb()x a+L sm+L train(7)L react()x a=T react⋅v()x a(8)L eb()x a=v2()x a2a eb(9)式中：S min()x a为列车在位置x a处的最小安全追踪171交通运输系统工程与信息2020年12月距离；L react()x a为列车在反应时间T react内的走行距离；L eb()x a为列车紧急制动距离；a eb为紧急制动率；L sm为安全余量；L train为列车车长.为简化计算，假设T react、a eb、L sm和L train为固定值.根据列车在通过最小安全追踪距离期间是否停靠站台，将最小追踪间隔的计算分为区间和车站两种模式.图1为区间最小追踪间隔示意，相邻两车以最小安全追踪距离为间隔分别运行至x a和x a+S min()x a处，其中，x a+S min()x a<x S，x S为站台位置.相邻列车沿图中实线所示的运行轨迹连续经过x a的最小时间间隔H min()x a为H min()x a=∫x a x a+S min()x a1v()x⋅d x(10)图1区间最小追踪间隔Fig.1Minimum headway of interstation地铁车站一般不设配线，列车到达车站后需在正线上完成停站作业，且同一时间只允许1列列车进行停站作业.车站最小追踪间隔的计算不考虑存在多配线条件下的列车到达、到通、出发、发通等间隔，如图2所示.相邻两车分别运行至x a和x a+S min()x a处，其中，x a+S min()x a≥x S，列车通过S min()x a期间需经历减速进站、停车和加速出站过程.因此，相邻列车连续经过x a的最小时间间隔H min()x a为H min()x a=T in()x a+T dwell+T out()x a(11)T in()x a=∫x a x S1v()x⋅d x(12)T out()x a=∫x S x a+S min()x a1v()x⋅d x(13)式中：T in()x a、T out()x a分别为列车进站、出站运行时间；T dwell为停站时间.图2车站区域最小追踪间隔Fig.2Minimum headway of station area对于任意位置x a，H min()x a取决于S min()x a和通过S min()x a的列车运行速度.由式(7)~式(9)可知，S min()x a取决于v()x a和x a以前的列车操纵策略有关，通过S min()x a的列车运行速度和x a以后的列车操纵策略有关.因此，操纵策略决定了列车追踪运行过程中可实现的最小追踪间隔.全线列车执行相同U在线路上追踪运行，可实现的最小追踪间隔H()U为列车连续通过最受限制点的最小时间间隔，即H()U=maxx∈[]0,XH min()x(14) 2列车运行性能优化2.1问题描述列车站间运行时间由运行图预先给定，理论上存在无数种操纵策略保证列车安全、正点运行.面向高峰小时地铁列车密集追踪的运行场景，将其中可实现追踪间隔最小且牵引能耗最小的操纵策略作为最优操纵策略.Takeuchi[5]对最小追踪间隔计算进行了灵敏度分析，结果表明，增大列车牵引/制动加速度可有效压缩最小追踪间隔.Liu[8]利用极大值原理推导出连续控制条件下的列车节能驾驶最优控制应包括最大牵引、巡航、惰行和最大制动.上述文献表明，列车最大牵引和最大制动是实现节能驾驶和高效追踪的必要控制条件.因此，本文选择最大牵引、巡航、惰行和最大制动组成列车最优操纵策略的4种控制变量，目标问题转化为寻找这些控制变量之间的组合顺序及其对应的转换点.本文假设列车装备自动驾驶系统，能够在线172第20卷第6期城市轨道交通列车追踪间隔与牵引能耗优化路任意位置实现巡航控制.因现有地铁系统并未完全装备再生制动设备，为不失一般性，再生制动能不在本文考虑范围内.2.2优化模型将U 作为决策变量，将min E ()U 和min H ()U 作为优化目标，设置列车安全、正点运行约束条件，构建双目标优化模型为ìíîïïïïmin E ()U min H ()U s.t.T ()U =T setv ()x <V limit ()x ,x ∈[0,X ](15)式中：V limit ()x 为线路限速条件；T set 为指定站间运行时间.2.3模型求解构建双目标优化模型旨在降低列车牵引能耗的同时压缩列车最小追踪间隔，期望得到一组准确的Pareto 最优解以体现两目标之间的均衡，采用基于动态规划方法的精确算法求解式(15).由于最小追踪间隔指标在动态规划逐段递推过程中并不严格单调，不具备动态规划方法所需的“无后效性”，引入ε-约束法将其转化为约束条件，并将式(15)转化为单目标优化模型，并以Δε为间隔由小到大调整ε参数值后多次利用动态规划方法进行求解.ìíîïïïïmin E ()U s.t .H ()U ≤εT ()U =T setv ()x <V limit ()x ,x ∈[0,X ](16)对式(16)的状态空间进行离散，将列车运行过程在空间域上划分K 个阶段，保证各阶段线路条件不变且各阶段长度不大于Δx .在k 阶段开始位置x k 处，按间隔Δv 划分k 阶段的开始状态集{}s k,i|sk,i=()x k ,i ⋅Δv ,0≤i ⋅Δv ≤Vˉ()x ,i ∈Ν，其中，Vˉ()x 为列车在最速策略下运行至x 处的速度.如图3所示，实线为列车最速运行轨迹，圆点为各阶段开始状态，亦为上一阶段的结束状态.列车运行过程的起点和终点分别记为s 1,0和s K +1,0.动态规划始终运行在基于上述离散划分规则的状态空间内，其优化结果一定能够满足安全运行约束v ()x <V limit ()x.图3模型阶段状态划分示意图Fig.3Stages and vertices of optimization model列车在k 阶段运行时固定采用最大牵引、巡航、惰行和最大制动中的一种工况，将其简记为u k ，列车在相邻阶段的最优控制工况切换还应满足图4所示的接续条件.图4列车控制接续约束Fig.4Train regime switching constraint列车在k 阶段s k,i 处施加u k 后运行至k +1阶段()x k +1,v ′处，()x k +1,v ′有可能不属于已划分的开始状态集，需要对其进行近似处理并修正至s k +1,j 处，其中j =round ()v ′Δv ，将上述状态转移过程简写为s k +1,j =F ()s k,i ,u k (17)列车在s k,i 处施加u k 阶段牵引能耗e ()s k,i ,u k 和运行时间t ()s k,i ,u k 计算为e ()s k,i ,u k =∫x kx k +1max []0,u ()x ⋅d x(18)t ()s k,i ,u k =∫xkx k +11v ()x ⋅d x(19)为获取式(16)在ε参数下的最优操纵策略U *ε，采用后向动态规划的求解方式，从s K +1,0开始逐段向前推进寻找列车在各阶段的最优控制决策173交通运输系统工程与信息2020年12月直至s 1,0.对于阶段k ，将列车从s k,i 运行至s K +1,0的操纵子策略记为U ()s k,i ，其对应的牵引能耗和运行时间分别为E []U () s k,i 和T []U ()s k,i ，下面建立评价U ()s k,i 最优性的过程指标函数.将列车从s 1,0运行至s k,i 处的操纵子策略记为U ()s k,i ，其对应的牵引能耗及运行时间分别记为E []U () s k,i 和T []U ()s k,i .因此，式(16)的优化目标min E ()U 可拆解为min E []U () s k,i 和min E []U ()s k,i .在后向动态规划求解过程中，E []U ()s k,i 无法直接求解得到.列车牵引能耗与运行时间在节能操纵条件下呈反比关系，故min E []U ()s k,i 等效于max T []U () s k,i .因T []U () s k,i +T []U ()s k,i =T set ，故min E []U () s k,i 等效于min T []U ()s k,i .因此，评价U ()s k,i 最优性的过程指标函数为E []U *() s k,i =min u k{}e ()s k,i ,u k +E {}U *éëùû F ()s k,i ,u k (20)T []U *() s k,i =min u k{}t ()s k,i ,u k +T {}U *éëùûF ()s k,i ,u k (21)U ()s k,i 对应的最小追踪间隔H []U () s k,i 应满足ε约束，即H []U ()s k,i ≤ε(22)U ()s k,i 应满足准点运行约束为T min []U () s k,i ≤T []U () s k,i ≤T max []U () s k,i (23)式中：T min []U () s k,i 、T max []U ()s k,i 分别为列车从s k,i 到s K +1,0的最短、最长运行时间.特别的，列车在s 1,0和s K +1,0的准点约束条件为T []U () s 1,0=T min []U () s 1,0=T max []U ()s 1,0=T set (24)T []U () s K +1,0=T min []U () s K +1,0=T max []U ()s K +1,0=0(25)T min []U () s k,i 和T max []U ()s k,i 的计算方法如图5所示.点划线为列车从s k,i 到s K +1,0的最速运行轨迹，其对应的运行时间为T min []U ()s k,i .虚线为列车从s 1,0到s k,i 的最速运行轨迹，其对应的运行时间为T min []U () s k,i ，则T max []U () s k,i =T set -T min []U () s k,i .图5T min []U () s k,i 和T max []U ()s k,i 计算原理Fig.5Calculation principle of T min []U () s k,i and T max []U ()s k,i由于存在2个冲突的过程指标，即式(20)和式(21)，列车从s k,i 处运行至s K +1,0的最优操纵子策略应为一组Pareto 解，记为{}U *w () s k,i ，其中，U *w ()s k,i 为s k,i 处的第w 个最优子策略.动态规划方法从s K +1,0开始逐段向前推进，求解各阶段各状态点的Pareto 最优子操纵策略集直至起点s 1,0.由于准点运行式(24)的限制，s 1,0处的最优操纵策略有且只有一个，为式(16)在ε参数下的最优解，即U *ε=U *1()s 1，0.式(16)的具体求解步骤描述如下：Step 1读入线路数据及列车参数并计算V ˉ()x ，设定模型参数Δx 、Δv 和ε，划分离散状态集合{}s k,i ，设定边界条件E []U *1()s K +1,0=0，T []U *1()s K +1,0=0，令k =K ，完成模型求解初始化.Step 2对于k 阶段所有的s k,i ，遍历4种最优控制工况u k ，获取k +1阶段F ()s k,i ,u k 处的最优子策略集{}U *w éëùûF ()s k,i ,u k .在满足式(22)和式(23)的条件下，将u k 加入{}U *w éëùûF ()s k,i ,u k 生成s k,i 处的可行子策略集{}U ()s k,i .Step 3根据过程指标式(20)和式(21)，对{}U ()s k,i 进行Pareto 占优操作，得到s k,i 的Pareto 最优子策略集{}U *w () s k,i .若k ≠1，则令k =k -1，跳转174第20卷第6期城市轨道交通列车追踪间隔与牵引能耗优化至Step 2；若k =1，输出式(16)的最优解U *ε.3算例分析基于C++开发优化程序，以北京地铁亦庄线为算例，选择文献[2]中线路数据和车辆参数，在此基础上进一步设定，a eb =-1.0m/s 2、L sm =30m 、L train =90m 、T react =0.5s 和T dwell =30s .根据“计算时间—优化效益可接受原则”设置模型参数如下：Δx =10m 、Δv =1km/h 和Δε=0.1，模型参数决定模型的求解精度和求解效率，间隔越小，模型求解精度越高且求解效率越低.以第7区间“万源街—荣京东”为例：首先，令ε=0，不断迭代计算ε=ε+Δε后的式(16)直至获取第1个有效解，该解即为能够实现理论最小追踪间隔的操纵策略U *headway ，令εmin =H ()U *headway ；其次，将ε设置为一个较大的正值后对式(16)进行求解，所得解即为只考虑节能目标的最优操纵策略U *energy ，令εmax =H ()U *energy ；最后，以Δε为间隔在[]εmin ,εmax 内由小到大调整ε参数值，多次求解式(16)，获得一组Pareto 最优解，如图6所示.图6列车最优操纵Pareto 解Fig.6Pareto solution of optimum driving strategy选取ε=70条件下的最优操纵策略U *70，对比U*energy、U *70和U*headway这3种策略下的v ()x 及其对应的H min ()x ，如图7所示.U*energy条件下列车从起点开始最大牵引至A 1，然后以惰行和巡航的组合方式运行至B 1，最后施加最大制动至终点；列车从C 1开始进入车站追踪模式，H min ()x 发生跃变；U *energy的性能指标分别为E ()U*energy=9.8kW ⋅h ，H ()U *energy =77.2s .U *70条件下列车最大牵引至A 2，惰行并巡航至B -2后开始施加第1次最大制动至C 2，然后惰行至B 2后施加第2次最大制动至终点；由于提前制动，列车在C 2处才进入车站追踪模式，其性能指标为E ()U *70=10.7kW∙h，H ()U *70=70.0s .同U *70一样，U *headway 在B -3处提前制动，采用两次制动的进站模式分别经过了A 3-B -3-C 3-B 3点，列车自C 3开始进入车站追踪模式，性能指标为E ()U *headway =13.0kW∙h ，H ()U *headway =67.1s .相较于U *energy，U *70和U *headway 在进站过程中执行两次制动策略，压缩了列车最小追踪间隔，为弥补两次制动进站过程中额外消耗的运行时间，U *70和U *headway 提升了进站以前的运行速度，额外增加了牵引能耗.图7U *energy 、U *70和U *headway 的列车运行轨迹及最小时间间隔Fig.7Train trajectory and minimum time separation ofU *energy ,U *70and U *headway列车进站操纵策略变化本质上影响的是车站追踪模式下T in ()x 和T out ()x ，从而改变H min ()x .U *energy 、U *70和U *headway 策略下的T in ()x 和T out ()x 随v ()x 的变化趋势如图8所示.上述操纵策略在车站追踪模式下的运行过程均包含1次惰行和1次最大制动.随着列车向终点运行，3种操纵策略对应的175交通运输系统工程与信息2020年12月T out()x从0开始递增并分别在列车进站的制动初始点B1、B2和B3处达到极大值，T in()x呈递减趋势并在终点处减为0.由于不同操纵策略下T in()x和T out()x存在差异，U*energy条件下的H min()x呈先增后减趋势，在B1处达到极大值；U*70和U*headway条件下，H min()x分别在C2和C3处达到极大值.图8车站追踪模式下列车运行轨迹及相关时间间隔Fig.8Train trajectory and related time separation under station tracking mode 计算全线13个站间的最优操纵Pareto解，选取各站间的U*headway与U*energy进行比较，如表1所示.与U*energy相比，列车在U*headway下的全线通行能力提升了17.0%，牵引总能耗增加了19.3%.各站间的U*headway均采用两次制动的进站策略且第2次制动的进站初速度都为23.0km/h，各站间U*headway的最小追踪间隔均达到约67s的极小值.各站间的U*energy都采用一次制动进站策略，其最小追踪间隔的大小与列车进站制动初速度正相关.不失一般性，图9以第2区间和第4区间为例，描绘U*headway和U*energy策略下v(x)和H min(x)来进一步印证上述结论.表1列车操纵策略优化结果Table1Optimaziton results of drivingstrategy176第20卷第6期城市轨道交通列车追踪间隔与牵引能耗优化图9U *headway 和U *energy 的列车运行轨迹及最小时间间隔Fig.9Train trajectory and minimu time separation of U *headway 和U *energy4结论列车操纵策略决定了列车在线路上运行的牵引能耗和可实现的最小追踪间隔.列车采用两次制动的进站策略可以有效压缩最小追踪间隔，但需要消耗更多的牵引能耗.列车最小追踪间隔和进站制动初速度正相关且存在极小值.相较于只考虑节能目标的最优操纵策略，追踪间隔压缩后的最优操纵策略最多可提升17.0%的线路通过能力，同时也增加了19.3%的牵引能耗.运营商可权衡地铁高峰小时服务水平和运营成本的实际需求，利用本文方法获得列车追踪运行最优操纵策略，具有一定实际意义.随着再生制动设备在城市轨道交通系统的应用普及，考虑再生制动能利用的列车运行能耗和追踪间隔多目标优化问题有待进一步讨论.参考文献：[1]WANG Y,DE SCHUTTER B,VAN DEN BOOM T J J,et al.Optimal trajectory planning for trains under fixed and moving signaling systems using mixed integer linear programming[J].Control Engineering Practice,2014,22:44-56.[2]YE H,LIU R.A multiphase optimal control method for multi-train control and scheduling on railway lines[J].Transportation Research Part B:Methodological,2016,93:377-393.[3]柏赟,于昭,贾文峥,等.考虑追踪安全的地铁快慢车协同操纵节能优化[J].交通运输系统工程与信息,2019,19(3):126-133.[BAI Y,YU Z,JIA W Z,et al.Cooperative control of express/local metro trains for energy saving considering the safe headway[J].Journal of Transportation Systems Engineering and Information Technology,2019,19(3):126-133.][4]WANG P,GOVERDE R M P.Multi-train trajectory optimization for energy-efficient timetabling[J].European Journal of Operational Research,2019,272(2):621-635.[5]TAKEUCHI H,GOODMAN C J,SONE S.Moving block signalling dynamics:Performance measures and re-starting queued electric trains[J].IEE Proceedings-electric Power Applications,2003,150(4):483-492.[6]陈荣武,诸昌钤,刘莉.CBTC 系统列车追踪间隔计算及优化[J].西南交通大学学报,2011,46(4):579-585.[CHEN R W,ZHU C Q,LIU L.Calculation and optimization of train headway in CBTC system[J].Journal of Southwest Jiaotong University,2011,46(4):579-585.]下转第204页177交通运输系统工程与信息2020年12月Influences of traffic flow characteristics on accidentseverity on secondary roads[J].China Journal ofHighway and Transport,2020,33(2):135-145.] [9]吴琴,施欣,陶学宗.海事事故严重性影响因素及影响程度识别[J].交通运输系统工程与信息,2019,19(1):189-195.[WU Q,SHI X,TAO X Z.Identifying the factors and their impact levels on severity of maritime traffic accidents[J].Journal of Transportation Systems Engineering and Information Technology,2019,19(1): 189-195.]上接第177页[7]NAKAMURA H.Analysis of minimum train headway ona moving block system by genetic algorithm[J].WITTransactions on the Built Environment,1998,37:1013-1022.[8]LIU R R,GOLOVITCHER I M.Energy-efficientoperation of rail vehicles[J].Transportation ResearchPart A:Policy and Practice,2003,37(10):917-932.上接第183页[2]李军,邓红平.基于公交IC卡数据的乘客出行分类研究[J].重庆交通大学学报(自然科学版),2016,35(6):109-114.[LI J,DENG H P.Classification ofpassenger's travel behavior based on IC card data[J].Journal of Chongqing Jiaotong University(NaturalSciences),2016,35(6):109-114.][3]何兆成,余畅,许敏行.考虑出行模式和周期性的公交出行特征分析[J].交通运输系统工程与信息,2016,16(6):135-141.[HE Z C,YU C,XU M X.Analyzingmethods of residents'travel characteristics consideringtravel patterns and periodicity[J].Journal ofTransportation Systems Engineering and InformationTechnology,2016,16(6):135-141.][4]GOULET-LANGLOIS G,KOUTSOPOULOS H N,ZHAO Z,et al.Measuring regularity of individual travelpatterns[J].IEEE Transactions on IntelligentTransportation Systems,2017,19(5):1583-1592. [5]翁剑成,王昌,王月玥,等.基于个体出行数据的公共交通出行链提取方法[J].交通运输系统工程与信息,2017,17(3):67-73.[WENG J C,WANG C,WANG YY,et al.Extraction method of public transit trip chainsbased on the individual riders'data[J].Journal ofTransportation Systems Engineering and InformationTechnology,2017,17(3):67-73.][6]ZOU Q,YAO X,ZHAO P,et al.Detecting home locationand trip purposes for cardholders by mining smart cardtransaction data in Beijing subway[J].Transportation,2018,45(3):919-944.[7]PEI J,HAN J,MORTAZAVI-ASL B,et al.Miningsequential patterns by pattern-growth:The prefixspanapproach[J].IEEE Transactions on Knowledge and DataEngineering,2004,16(11):1424-1440.204。

第９卷㊀第２期２０２４年３月气体物理ＰＨＹＳＩＣＳＯＦＧＡＳＥＳＶｏｌ.９㊀Ｎｏ.２Ｍａｒ.２０２４㊀㊀ＤＯＩ:１０.１９５２７/ｊ.ｃｎｋｉ.２０９６￣１６４２.１０８９可压缩两气体流动的简化神经网络模型刘子岩ꎬ㊀许㊀亮ꎬ㊀刘耀峰(中国航天空气动力技术研究院ꎬ北京１０００７４)ＡＳｉｍｐｌｉｆｉｅｄＮｅｕｒａｌＮｅｔｗｏｒｋＭｏｄｅｌｆｏｒＣｏｍｐｒｅｓｓｉｂｌｅＴｗｏ￣ＧａｓＦｌｏｗｓＬＩＵＺｉｙａｎꎬ㊀ＸＵＬｉａｎｇꎬ㊀ＬＩＵＹａｏｆｅｎｇ(ＣｈｉｎａＡｃａｄｅｍｙｏｆＡｅｒｏｓｐａｃｅＡｅｒｏｄｙｎａｍｉｃｓꎬＢｅｉｊｉｎｇ１０００７４ꎬＣｈｉｎａ)摘㊀要:实用的虚拟流体方法(ｐｒａｃｔｉｃａｌｇｈｏｓｔｆｌｕｉｄｍｅｔｈｏｄꎬＰＧＦＭ)利用Ｒｉｅｍａｎｎ问题速度解对可压缩多介质流场界面条件进行建模ꎮ基于构造的嵌入物理约束的神经网络模型预测Ｒｉｅｍａｎｎ问题速度解的方式ꎬ给出一种两气体流动的神经网络模型简化方法ꎮ首先提出完全气体状态方程下神经网络模型输入特征采样范围从无界域到有界域的转换方法ꎬ改善模型预测不同初始条件下Ｒｉｅｍａｎｎ解的泛化性能ꎮ根据该转化方法ꎬ进一步提出一种结构更加简单的神经网络优化方法ꎬ将输入维度从５个减少到３个ꎬ有效提高神经网络的训练效果ꎮ将该神经网络代理模型应用于ＰＧＦＭ程序框架ꎬ通过典型的一维与二维两气体流动问题进行数值验证与对比分析ꎮ结果表明ꎬ简化的网络模型与已有研究的神经网络模型相比ꎬ能取得精度相近的计算结果ꎮ而在神经网络训练效率上ꎬ简化神经网络具有明显优势ꎮ同时因为简化神经网络采样维度少ꎬ方便尝试加密采样提高拟合精度ꎬ更具备发展潜力ꎮ关键词:可压缩多介质问题ꎻ虚拟流体方法ꎻ两气体Ｒｉｅｍａｎｎ问题ꎻ神经网络㊀㊀㊀收稿日期:２０２３￣０９￣２６ꎻ修回日期:２０２３￣１２￣１９基金项目:国家自然科学基金(１１８７２３５１)第一作者简介:刘子岩(１９９８ )㊀男ꎬ博士ꎬ主要研究方向为机器学习在可压缩多介质流问题中的应用ꎮＥ￣ｍａｉｌ:ｌｚｙ２１０＠ｂｕａａ.ｅｄｕ.ｃｎ通信作者简介:许亮(１９８２ )㊀男ꎬ研究员ꎬ主要研究方向为可压缩多介质流模拟方法㊁ＣＦＤ中机器学习建模方法等ꎮ中图分类号:Ｏ３５９㊀㊀文献标志码:ＡＡｂｓｔｒａｃｔ:Ｔｈｅｐｒａｃｔｉｃａｌｇｈｏｓｔｆｌｕｉｄｍｅｔｈｏｄ(ＰＧＦＭ)ｕｔｉｌｉｚｅｓｖｅｌｏｃｉｔｙｓｏｌｕｔｉｏｎｓｏｆＲｉｅｍａｎｎｐｒｏｂｌｅｍｓｔｏｍｏｄｅｌｔｈｅｉｎｔｅｒｆａｃｅｅｖｏｌｕｔｉｏｎｏｆｃｏｍｐｒｅｓｓｉｂｌｅｍｕｌｔｉ￣ｍａｔｅｒｉａｌｆｌｏｗｓ.Ｔｈｉｓｐａｐｅｒｐｒｅｓｅｎｔｅｄａｓｉｍｐｌｉｆｉｅｄｎｅｕｒａｌｎｅｔｗｏｒｋｍｏｄｅｌｆｏｒｔｗｏ￣ｇａｓｆｌｏｗｓｂｙｐｒｅｄｉｃｔｉｎｇｔｈｅｖｅｌｏｃｉｔｙｓｏｌｕｔｉｏｎｏｆＲｉｅｍａｎｎｐｒｏｂｌｅｍｂａｓｅｄｏｎｔｈｅｎｅｕｒａｌｎｅｔｗｏｒｋｍｏｄｅｌｅｍｂｅｄｄｅｄｗｉｔｈｐｈｙｓｉｃａｌｃｏｎｓｔｒａｉｎｔｓ.Ｆｉｒｓｔｌｙꎬａｍｅｔｈｏｄｆｏｒｃｏｎｖｅｒｔｉｎｇｔｈｅｓａｍｐｌｉｎｇｒａｎｇｅｏｆｔｈｅｎｅｕｒａｌｎｅｔｗｏｒｋｍｏｄｅｌｆｒｏｍｕｎｂｏｕｎｄｅｄｄｏｍａｉｎｔｏｂｏｕｎｄｅｄｄｏｍａｉｎｗａｓｐｒｏｐｏｓｅｄꎬｗｈｉｃｈｈｏｌｄｓｔｒｕｅｆｏｒｔｈｅｐｅｒｆｅｃｔｇａｓｅｑｕａｔｉｏｎｏｆｓｔａｔｅ.Ｉｔｃａｎｉｍｐｒｏｖｅｔｈｅｇｅｎｅｒａｌｉｚａｔｉｏｎｐｅｒｆｏｒｍａｎｃｅｏｆｔｈｅｍｏｄｅｌｕｎｄｅｒｄｉｆｆｅｒｅｎｔｉｎｉｔｉａｌｃｏｎｄｉｔｉｏｎｓ.Ｂａｓｅｄｏｎｔｈｉｓｔｒａｎｓｆｏｒｍａｔｉｏｎｍｅｔｈｏｄꎬａｓｉｍｐｌｅｒｎｅｕｒａｌｎｅｔｗｏｒｋｓｔｒｕｃｔｕｒｅｗａｓｆｕｒｔｈｅｒｐｒｏｐｏｓｅｄ.Ｔｈｅｔｒａｉｎｉｎｇｒｅｓｕｌｔｏｆｔｈｅｎｅｕｒａｌｎｅｔｗｏｒｋｃａｎｂｅｅｆｆｅｃｔｉｖｅｌｙｉｍｐｒｏｖｅｄｂｙｒｅｄｕｃｉｎｇｔｈｅｉｎｐｕｔｄｉｍｅｎｓｉｏｎｓｆｒｏｍ５ｔｏ３.ＴｈｅｎｅｕｒａｌｎｅｔｗｏｒｋｍｏｄｅｌｗａｓａｐｐｌｉｅｄｔｏｔｈｅＰＧＦＭ.Ｎｕｍｅｒｉｃａｌｖａｌｉｄａｔｉｏｎｏｆｔｈｅｎｅｕｒａｌｎｅｔｗｏｒｋｍｏｄｅｌｗａｓｃａｒｒｉｅｄｏｕｔｔｈｒｏｕｇｈｔｙｐｉｃａｌｏｎｅ￣ｄｉｍｅｎｓｉｏｎａｌａｎｄｔｗｏ￣ｄｉｍｅｎｓｉｏｎａｌｇａｓｆｌｏｗｐｒｏｂｌｅｍｓ.Ｔｈｅｒｅｓｕｌｔｓｓｈｏｗｔｈａｔｔｈｅｓｉｍｐｌｉｆｉｅｄｎｅｔｗｏｒｋｍｏｄｅｌｃａｎａｃｈｉｅｖｅｓｉｍｉｌａｒｃｏｍｐｕｔａｔｉｏｎａｌａｃｃｕｒａｃｙｃｏｍｐａｒｅｄｗｉｔｈｅｘｉｓｔｉｎｇｎｅｕｒａｌｎｅｔｗｏｒｋｍｏｄｅｌｓ.Ｉｎｔｅｒｍｓｏｆｔｒａｉｎｉｎｇｅｆｆｉｃｉｅｎｃｙｏｆｎｅｕｒａｌｎｅｔｗｏｒｋｓꎬｔｈｅｓｉｍｐｌｉｆｉｅｄｎｅｕｒａｌｎｅｔｗｏｒｋｈａｓｏｂｖｉｏｕｓａｄｖａｎｔａｇｅｓ.Ｍｏｒｅｏｖｅｒꎬｂｅｃａｕｓｅｔｈｅｓｉｍｐｌｉｆｉｅｄｎｅｕｒａｌｎｅｔｗｏｒｋｈａｓｆｅｗｅｒｓａｍｐｌｉｎｇｄｉｍｅｎｓｉｏｎｓꎬｉｔｉｓｃｏｎｖｅｎｉｅｎｔｔｏｔｒｙｄｅｎｓｅｒｓａｍｐｌｉｎｇｔｏｉｍｐｒｏｖｅｆｉｔｔｉｎｇａｃｃｕｒａｃｙａｎｄｓｕｃｈｍｅｔｈｏｄｈａｓｍｏｒｅｄｅｖｅｌｏｐｍｅｎｔｐｏｔｅｎｔｉａｌ.Ｋｅｙｗｏｒｄｓ:ｃｏｍｐｒｅｓｓｉｂｌｅｍｕｌｔｉ￣ｍａｔｅｒｉａｌｆｌｏｗｓꎻｇｈｏｓｔｆｌｕｉｄｍｅｔｈｏｄꎻｔｗｏ￣ｇａｓＲｉｅｍａｎｎｐｒｏｂｌｅｍꎻｎｅｕｒａｌｎｅｔｗｏｒｋ气体物理２０２４年㊀第９卷引㊀言数值模拟可压缩多介质流体动力学问题在航空航天㊁武器物理㊁生物医疗等前沿领域具有大量应用需求[１]ꎮ清晰准确地模拟界面运动及变形是这类问题关注的重点ꎮ处理物质界面的一类方案为锐利界面方法ꎬ即把界面视作零厚度界面ꎬ这种方案须时刻捕捉界面位置并定义非线性波与界面发生作用而产生的边界条件ꎮ捕捉界面位置常用的方法有界面追踪法[２]和水平集(ｌｅｖｅｌｓｅｔ)方法[３]等ꎬ定义边界条件的方法有切割网格法[４]和虚拟流体方法等ꎮ虚拟流体方法有着诸多版本ꎬ其中Ｌｉｕ等[５]提出的修正虚拟流体方法(ｍｏｄｉｆｉｅｄｇｈｏｓｔｆｌｕｉｄｍｅｔｈｏｄꎬＭＧＦＭ)通过求解界面附近多介质Ｒｉｅｍａｎｎ问题来预测界面状态ꎬ并以此定义虚拟流体状态ꎮ许亮等[６]提出了实用虚拟流体方法ꎬ借助多介质Ｒｉｅｍａｎｎ问题的速度解ꎬ并结合界面另一侧的流场状态可以准确定义虚拟流体状态ꎮ由于充分考虑波与物质界面的相互作用和介质性质对界面运动的影响ꎬ这两种方法已被证明适用于各种强间断气－气㊁气－液㊁流－固等挑战性多介质流动难题[７]ꎮ使用基于Ｒｉｅｍａｎｎ解的虚拟流体方法在定义虚拟流体状态时ꎬ须迭代求解非线性方程ꎮ多介质Ｒｉｅｍａｎｎ问题近似求解可以分成隐式和显式两类方法[８]ꎮ虚拟流体方法中常用的基于双激波结构的隐式近似Ｒｉｅｍａｎｎ解[５]也包含迭代过程ꎮＸｕ等[８]发展了２种不含迭代的显式近似Ｒｉｅｍａｎｎ解ꎬ计算效率高ꎬ但是只适用于刚性气体状态方程ꎮ传统方法求解多介质Ｒｉｅｍａｎｎ问题应对不同的状态方程须进行特定的数学推导ꎮ当问题比较复杂时ꎬ通常会涉及繁琐的计算过程ꎮ比如ꎬ对于含热完全气体的多介质流动问题ꎬ在压力的每一步迭代求解过程中也包含了温度的迭代解[８]ꎮ如果考虑黏性㊁表面张力㊁相变㊁化学反应等影响因素ꎬ那么物质界面求解迭代过程可能更加复杂ꎬ这就对稳定性提出很高要求ꎮ因此ꎬ有必要尝试借助新兴的机器学习方法替代传统的Ｒｉｅｍａｎｎ解算器ꎬ对降低建模难度以及提高计算稳定性等方面进行探索ꎮ机器学习方法近年来被广泛应用于各种研究领域中[９]ꎬ尤其常见于求解偏微分方程的过程中替代复杂的计算ꎬ并取得了较好的效果[１０ꎬ１１]ꎮＦｅｎｇ等[１２ꎬ１３]使用神经网络构造了智能激波指示器ꎬ应用于格式设计ꎬ并给出了机器学习模型的合理性解释ꎮ同时ꎬ研究者们也尝试加入一些约束条件以使得结果更加符合流动物理ꎮ丘润荻等[１４]在相场法中引入含物理约束的神经网络ꎬ用于求解ＲＴ不稳定性问题ꎬ明显提高了计算效率ꎮＭａｇｉｅｒａ等[１５]针对Ｒｉｅｍａｎｎ问题ꎬ首次提出嵌入物理约束的神经网络模型ꎬ使预测结果尽可能满足Ｒａｎｋｉｎｅ￣Ｈｕｇｏｎｉｏｔ间断条件ꎮ基于该方法ꎬＷａｎｇ等[１６]针对跨临界流的Ｒｉｅｍａｎｎ问题ꎬ分析了具有不同输入与输出形式的神经网络架构ꎬ提出了兼顾计算效率及计算精度的网络构造方法ꎮＲｕｇｇｅｒｉ等[１７]将Ｒｉｅｍａｎｎ解的神经网络模型推广到真实气体ꎬ比精确Ｒｉｅｍａｎｎ解具有更高的计算效率ꎮ但是ꎬ由于关注的网络输出维度较多ꎬ这些文献所使用的神经网络结构仍较复杂ꎮ本人在文献[９]中提出了多介质Ｒｉｅｍａｎｎ问题速度解的神经网络模型ꎬ通过构建物理约束层的方式考虑流场间断关系ꎬ简化网络规模ꎬ成功应用于ＰＧＦＭ框架求解多介质流动问题ꎮ本文将进一步提出应用于两气体Ｒｉｅｍａｎｎ问题的神经网络标准化映射方法ꎬ构造针对输入特征的神经网络结构优化方式ꎬ能够提高网络的训练效果及应用价值ꎮ１㊀问题描述与求解方法１.１㊀流体控制方程与状态方程研究模型设定为可压缩无黏流动模型ꎬ其二维形式Ｅｕｌｅｒ守恒方程组为∂Ｕ∂ｔ＋∂Ｆ∂ｘ＋∂Ｇ∂ｙ＝０(１)式中ꎬＵ＝[ρꎬρｕꎬρｖꎬＥ]Ｔ代表守恒变量ꎬＦ＝[ρｕꎬρｕ２＋ｐꎬρｕｖꎬ(Ｅ＋ｐ)ｕ]ＴꎬＧ＝[ρｖꎬρｕｖꎬρｖ２＋ｐꎬ(Ｅ＋ｐ)ｖ]Ｔ均为守恒通量ꎮρꎬｐꎬｕ分别是流体的密度ꎬ压力以及速度ꎬＥ是流体总能量ꎬＥ＝ρｅ＋ρｕ２/２ꎮｅ是单位质量流体的内能ꎮ上述方程组的封闭还须补充一个状态方程ꎮ本文目前的工作只涉及由完全气体状态方程ｐ＝(γ－１)ρｅ(２)描述的可压缩两气体流动问题ꎮ式中ꎬγ表示气体比热比ꎮ１.２㊀虚拟流体方法在虚拟流体方法中ꎬ把物质界面看成内边界ꎮ一般每种介质取界面附近３~５个网格点的区域作４３第２期刘子岩ꎬ等:可压缩两气体流动的简化神经网络模型为虚拟流体区域ꎬ虚拟流体层数由计算采用的数值格式决定ꎮ该区域内的虚拟流体状态作为内边界条件来定义ꎮ这样可以将多介质界面问题转化成多个单介质问题来模拟ꎮ原则上ꎬ只要这些虚拟流体状态定义准确ꎬ任何高分辨率计算格式可以直接使用ꎬ不需要在界面附近特殊处理ꎮ本文使用水平集方法捕捉界面运动ꎬ使用ＰＧ￣ＦＭ[７]定义不同流体的界面条件ꎬＰＧＦＭ的原理如图１所示ꎮ对比ＭＧＦＭ[５]须使用３个界面自由度定义虚拟流体状态ꎬＰＧＦＭ仅借助１个界面速度就可以计算出虚拟流体速度ꎬ虚拟流体的压力和密度通过镜像对称得到ꎮ这种区别意味着当采用神经网络对流场界面条件建模时ꎬ在ＰＧＦＭ的框架下仅关注界面速度这一个物理量ꎬ神经网络的输入和输出更加简洁ꎬ并能因此取得更好的训练效果ꎮ图１㊀ＰＧＦＭ定义虚拟流体状态的示意图[７]Ｆｉｇ.１㊀ＩｌｌｕｓｔｒａｔｉｏｎｏｆｔｈｅＰＧＦＭｆｏｒｄｅｆｉｎｉｎｇｇｈｏｓｔｆｌｕｉｄｓｔａｔｅｓ[７]为抑制界面附近的非物理解ꎬ界面附近和虚拟流体的密度也可以通过某些校正技术[７]获得ꎮ１.３㊀多介质Ｒｉｅｍａｎｎ问题的速度解在多维问题中ꎬ虚拟流体方法首先须沿界面法向构造并求解多介质Ｒｉｅｍａｎｎ问题ꎮ一维形式的初始值间断Ｒｉｅｍａｎｎ问题为∂Ｕ∂ｔ＋∂Ｆ∂ｘ＝０ꎬＵ｜ｔ＝０＝ＵＬꎬｘ<ｘＩＵＲꎬｘ>ｘＩ{(３)其中ꎬＵ＝[ρꎬρｕꎬＥ]ＴꎬＦ＝[ρｕꎬρｕ２＋ｐꎬ(Ｅ＋ｐ)ｕ]ＴꎬＵＬ和ＵＲ是界面法向上由位于ｘＩ处的物质界面分开的两个常值状态ꎮ下标Ｉ㊁Ｌ和Ｒ分别表示界面㊁左侧和右侧ꎮ对多介质Ｒｉｅｍａｎｎ问题(３)ꎬ须求解下面的非线性系统ｕＬ－ｕＩ＝ｆＬ(ｐＩ)ｕＩ－ｕＲ＝ｆＲ(ｐＩ){(４)其中函数ｆ的具体形式为[６]ｆＨ(ｐ)＝ｆ(ｐＩꎬｐＨꎬρＨ)＝(ｐ－ｐＨ)/(ｐ－ｐＨ)ρＨρρ－ρＨꎬｐ>ｐＨʏｐｐＨ１ρｄρｄｐｄｐꎬｐɤｐＨìîíïïïïïï(５)式中ꎬＨ＝Ｌ或Ｒꎬρ和ｐ的关系由状态方程确定ꎮ由ｕＩ＝１２(ｆＲ－ｆＬ)＋１２(ｕＬ＋ｕＲ)(６)求出界面速度解ｕＩꎬ然后按照图１中ＰＧＦＭ方式确定虚拟流体状态Ｕ∗Ｒ和Ｕ∗Ｌꎮ这样ꎬ每个时间步在界面处只需求解如下两个带虚拟流体状态的单介质Ｒｉｅｍａｎｎ问题∂Ｕ∂ｔ＋∂Ｆ∂ｘ＝０ꎬＵ｜ｔ＝０＝ＵＬꎬｘ<ｘ０Ｕ∗Ｒꎬｘ>ｘ０{(７)∂Ｕ∂ｔ＋∂Ｆ∂ｘ＝０ꎬＵｔ＝０＝Ｕ∗Ｌꎬｘ<ｘ０ＵＲꎬｘ>ｘ０{(８)这里上标∗表示虚拟流体状态ꎮ问题(７)的求解区域是从最左边的网格点一直跨过物质界面进入右边介质区的虚拟网格点ꎮ问题(８)的求解区域则是从左边的虚拟网格点跨过物质界面一直到最右边的网格点ꎮ２㊀神经网络模型建立及优化２.１㊀嵌入物理约束的神经网络模型神经网络的输入和输出须根据具体问题确定ꎮＰＧＦＭ定义虚拟流体状态只需利用多介质Ｒｉｅｍａｎｎ问题界面速度解ｕＩꎬ而根据界面速度ｕＩ的求解流程(４)~(６)ꎬ神经网络的输入特征可选择为{ｐＬꎬρＬꎬｐＲꎬρＲꎬΔｕ}共５个量ꎬ神经网络输出特征为函数(５)ｆＬ和ｆＲ的逼近值φＬ和φＲꎮ最终根据神经网络预测结果和式(６)计算得到ｕＩꎮ为了使神经网络学习结果更符合流动物理要求ꎬ根据流场间断关系构建了物理约束条件φＬ＋φＲ＋Δｕ＝０(９)通过构建物理约束层的方式将式(９)嵌入神经网络ꎬ具体操作方法为:将神经网络两个输出φＬ和φＲ改写成φＬ和－Δｕ－φＬꎬ神经网络模型的输出层只需要一个输出量φＬꎮ但神经网络训练样本的标签仍为两个量ꎬ由φＬ和－Δｕ－φＬ物理约束层参与训练ꎮ这种方式能严格满足式(９)的流场间断关５３气体物理２０２４年㊀第９卷系ꎮ考虑物理约束的神经网络原理如图２所示ꎮ图２㊀考虑物理约束的人工神经网络示意图Ｆｉｇ.２㊀Ｉｌｌｕｓｔｒａｔｉｏｎｏｆｔｈｅａｒｔｉｆｉｃｉａｌｎｅｕｒａｌｎｅｔｗｏｒｋｗｉｔｈｐｈｙｓｉｃａｌｃｏｎｓｔｒａｉｎｔｓ本部分工作已在文献[９]中完成ꎬ具体实现方法和神经网络的基本原理参考文献[９]ꎮ２.２㊀神经网络标准化映射及网络结构简化神经网络代理模型的使用范围应尽量在训练样本采样范围内以保证模型的预测效果ꎮ注意到ＰＧＦＭ的神经网络模型的５个输入量都是无界的ꎮ本质上ꎬ流体的压力和密度取值范围为(０ꎬ＋ɕ)ꎬ界面两侧流体的速度差Dｕ取值范围为(－ɕꎬ＋ɕ)ꎮ为了提高神经网络模型的训练效果ꎬ通过预处理的方法将采样范围从无界域映射到有界域ꎮ基于Ｒｉ￣ｅｍａｎｎ速度解的求解方法可以获得标准化映射方法依赖紧密的一个性质[１８]ꎮ对于完全气体状态方程(２)ꎬ如果界面两侧流体的压力㊁密度经过以下变换ｐ－＝ｐｐ０ꎬρ－＝ρρ０那么界面速度ｕＩ的表达式将被改写成ｕＩ＝１２[ｆ(ｐ－Ｉꎬｐ－Ｒꎬρ－Ｒ)－ｆ(ｐ－Ｉꎬｐ－Ｌꎬρ－Ｌ)]ｐ０ρ０＋１２(ｕＬ＋ｕＲ)(１０)其中ꎬｐ０ꎬρ０为任意大于零的参考压力和参考密度ꎮ为简化网络结构的输入输出关系ꎬ提高神经网络训练效果ꎬ降低采样难度ꎬ考虑将参考压力和参考密度设置为左右两侧压力㊁密度的较大值ꎬ即ｐ０＝ｍａｘ(ｐＬꎬｐＲ)ꎬρ０＝ｍａｘ(ρＬꎬρＲ)显然有０<ｐｍｉｎｐ０ɤ１ꎬｐｍａｘｐ０＝１ꎻ０<ρｍｉｎρ０ɤ１ꎬρｍａｘρ０＝１(１１)经过这样处理ꎬ５个输入特征中的２个特征确定为１ꎬ采样及训练时只须考虑其他３个不为１的输入特征ꎬ达到了简化网络结构的目的ꎮ值得注意的是ꎬ本方法研究对象为两气体问题ꎬ界面两侧流体的不同导致不能通过训练一个神经网络解决该问题ꎮ根据流体左右压力㊁密度的大小关系ꎬ共计有４种情况ꎬ须对应训练４个神经网络代理模型并综合使用ꎬ见表１ꎮ简化的神经网络标准化原理如图３ꎬ该图以ｐＬ<ｐＲꎬρＬȡρＲ的情况为例ꎬ神经网络选用２号神经网络ꎮ表１㊀各神经网络代理模型Ｔａｂｌｅ１㊀ＶａｒｉｏｕｓｎｅｕｒａｌｎｅｔｗｏｒｋｍｏｄｅｌｓＮｏ.ｒｅｌａｔｉｏｎｏｆｐｒｅｌａｔｉｏｎｏｆρｎｅｕｒａｌｎｅｔｗｏｒｋｉｎｐｕｔｓ１ｐＬ<ｐＲρＬ<ρＲ{ｐ－Ｌꎬρ－ＬꎬΔｕ－}２ｐＬ<ｐＲρＬȡρＲ{ｐ－Ｌꎬρ－ＲꎬΔｕ－}３ｐＬȡｐＲρＬȡρＲ{ｐ－Ｒꎬρ－ＲꎬΔｕ－}４ｐＬȡｐＲρＬ<ρＲ{ｐ－Ｒꎬρ－ＬꎬΔｕ－}图３㊀简化神经网络标准化原理示意图Ｆｉｇ.３㊀Ｉｌｌｕｓｔｒａｔｉｏｎｏｆｔｈｅｎｅｕｒａｌｎｅｔｗｏｒｋａｆｔｅｒｓｉｍｐｌｉｆｉｃａｔｉｏｎ因为每次计算都只会用到其中一个神经网络代理模型ꎬ相比于之前工作的神经网络模型[９]ꎬ简化后的神经网络模型降低了建模复杂度ꎬ同时也减小了计算量ꎬ可以提高计算效率ꎮ另外ꎬ根据前文介绍的性质ꎬ神经网络的输入经过了标准化映射ꎬ输出值也会变成标准化后的结果φ－Ｌ和φ－Ｒꎬ其中φ－Ｒ＝－Δｕ－－φ－ＬꎬΔｕ－＝Δｕ/ｐ０ρ０ꎬ为标准化后的速度差ꎮ须经过一个逆向映射得到φＬ６３第２期刘子岩ꎬ等:可压缩两气体流动的简化神经网络模型和φＲꎬ表达式为φＬ＝ｕ０φ－ＬꎬφＲ＝ｕ０φ－Ｒ(１２)ｕＩ即可根据式(１３)求出ｕＩ＝１２(φＲ－φＬ)＋１２(ｕＬ＋ｕＲ)(１３)２.３㊀训练样本生成神经网络模型的训练集通过一维多介质Ｒｉｅｍａｎｎ精确解算器生成ꎬ即求解非线性系统(４)ꎮ一般根据状态方程的具体表达形式ꎬ通过迭代方法获得界面压力解ｐＩꎬ再根据函数ｆＨ(ｐ)表达形式(５)获得神经网络的样本标签ｆＬ和ｆＲꎮ对于本文关注的完全气体状态方程(２)ꎬ其压力迭代方式和求解流程可以参考文献[１９]ꎮ采样范围根据２.２节的标准化映射方式而定ꎬ由式(１１)ꎬ界面两侧流体压力的较大值映射为１ꎬ较小值映射为一个大于０㊁小于等于１的值ꎻ界面两侧流体密度的较大值映射为１ꎬ较小值映射为一个大于０㊁小于等于１的值ꎮ所以只需对压力㊁密度的较小值进行采样ꎮ本文采用均匀采样的方式ꎬ界面两侧流体的压力㊁密度的较小值在(０ꎬ１]的范围内各取１１个点ꎮ标准化后的速度差Δｕ－与ｐ０和ρ０的选取相关ꎬ无法映射到一个具体的有界范围内ꎬ因此需要选定一个较为合理的范围ꎮ在实际工况下ꎬ压力ｐ０通常在１０５Ｐａ的量级上ꎬ密度ρ０在１ｋｇ/ｍ３的量级上ꎮ如果速度差的采样范围取在[－１０ꎬ１０]ꎬ即可覆盖[－３０００ｍ/ｓꎬ３０００ｍ/ｓ]的范围ꎬ这能满足绝大部分实际工况的需求ꎮ所以速度差在[－１０ꎬ１０]的范围内取２１个点ꎮ对样本数据的初始条件精确求解Ｒｉｅｍａｎｎ问题ꎬ并根据压力正项条件(Δｕ)ｃｒｉｔʉ２ｃＬγＬ－１＋２ｃＲγＲ－１>Δｕ(１４)去除其中的非物理解ꎬ将结果依次输出ꎬ即可得到训练神经网络需要的样本集ꎮ其中ꎬ(Δｕ)ｃｒｉｔ表示临界的速度差ꎬｃＬ和ｃＲ表示左右声速ꎮ减少神经网络输入维度带来的最直接的好处是大幅压缩了样本量ꎮ简化前的网络结构共有５个输入量ꎬ在忽略非物理解数量的情况下共计１１ˑ１１ˑ１１ˑ１１ˑ２１＝３０７４６１个样本ꎮ而３个输入量的神经网络要达到同样的采样密度ꎬ只需要使用１１ˑ１１ˑ２１＝２５４１个样本ꎮ同时ꎬ简化的神经网络拟合难度更低ꎬ这将大大减少神经网络的训练耗时ꎬ提升该方法的应用价值ꎮ神经网络的训练优势及计算效果将在２.４节算例应用部分进行展示ꎮ２.４㊀神经网络模型的应用方式４个神经网络完成训练后将其相关参数输出ꎬ以数学表达式的形式编制成程序模块嵌入到ＰＧＦＭ计算框架中[９]ꎬ即可直接应用ꎮ对于某个具体的两气体Ｒｉｅｍａｎｎ问题ꎬ首先判断界面左侧的压力㊁密度大小关系ꎬ标准化后应用对应神经网络代理模型ꎬ计算得到φＬ和φＲꎮ而后求解出界面速度ｕＩꎬ即可通过ＰＧＦＭ定义虚拟流体状态ꎮ类似于反射边界条件ꎬ对于一维问题ꎬ定义界面左侧流场的虚拟流体状态为ｐ∗＝ｐＲ㊀ρ∗＝ρＲ㊀ｕ∗＝２ｕＩ－ｕＲ(１５)定义界面右侧流场的虚拟流体状态为ｐ∗＝ｐＬ㊀ρ∗＝ρＬ㊀ｕ∗＝２ｕＩ－ｕＬ(１６)上标∗表示虚拟节点处的状态量ꎮ对于多维问题ꎬ需要根据ＰＧＦＭ的原理ꎬ将预测的界面法向速度与已知的切向速度合成ꎬ构造虚拟流体速度分量ꎬ保持密度和压力的定义不变ꎮ定义界面左侧流场的虚拟流体状态为ｐ∗Ｌ＝ｐＲꎬρ∗Ｌ＝ρＲꎬＶ∗Ｌ＝ＶＲ＋２(ｕｎꎬＩ－ｕｎꎬＲ)ｎ(１７)定义界面右侧流场的虚拟流体状态为ｐ∗Ｒ＝ｐＬꎬρ∗Ｒ＝ρＬꎬＶ∗Ｒ＝ＶＬ＋２(ｕｎꎬＩ－ｕｎꎬＬ)ｎ(１８)这里的下标ｎ表示界面法线方向ꎮ３㊀算例验证３.１㊀训练误差收敛效果本算例通过对比神经网络训练过程中误差的收敛速度和精度体现简化的神经网络结构的优势ꎮ假设左侧气体为空气(γＬ＝１.４)ꎬ右侧气体为氦气(γＲ＝１.６６７)ꎮ通过对比实验发现ꎬ神经网络隐含层层数为２ꎬ即可取得不错的拟合效果ꎮ为表述方便ꎬ本文使用Ｍ１￣Ｍ２￣Ｍ３￣Ｍ４的形式表示神经网络模型的规模ꎬ例如５￣１０￣１０￣２表示的是该神经网络模型输入量为５ꎬ输出量为２ꎬ含２个隐含层ꎬ每个隐含层各包含１０个节点ꎮ本文中所有训练均采用ＰｙＴｏｒｃｈ神经网络框架ꎬ优化过程使用ａｄａｍ优化器ꎮ７３气体物理２０２４年㊀第９卷图４分别给出了隐含层包含５个节点的神经网络模型和包含４０个节点的神经网络模型的两种网络结构训练效果对比ꎮ因为两种网络结构的训练集规模不同且相差较大ꎬ所以收敛对比图的横轴选为训练批量(ｂａｔｃｈ)数而不是训练周期(ｅｐｏｃｈ)数ꎻ纵轴ｌｏｓｓ为神经网络预测值与标签值的均方误差ꎮ具体计算式为ＬＭＳＥ＝１Ｎð ｙ－Ｆθ(ｘ) ２(１９)其中ꎬθ＝{ｗｌꎬｂｌ}Ｋｌ＝１为神经网络各层权重和偏置参数ꎬＦθ(ｘ)为神经网络预测结果ꎬｙ为标签值ꎬＮ表示每个训练批量的样本个数ꎮ训练过程中学习率设定为０.０００２ꎮ从图中结果可以反映出ꎬ３个输入量的神经网络模型在收敛速度上存在明显优势ꎮ５个输入量的神经网络模型也能达到同样的收敛精度ꎬ但是还需要花费相当多的训练步数ꎬ并且误差曲线依然存在明显的振荡ꎮ图中仅分别展示２００００和５０００个批量下的训练过程对比ꎮ(ａ)Ｎｅｕｒａｌｎｅｔｗｏｒｋｗｉｔｈ５ｎｏｄｅｓｉｎｔｈｅｈｉｄｄｅｎｌａｙｅｒ(ｂ)Ｎｅｕｒａｌｎｅｔｗｏｒｋｗｉｔｈ４０ｎｏｄｅｓｉｎｔｈｅｈｉｄｄｅｎｌａｙｅｒ图４㊀神经网络训练误差收敛图Ｆｉｇ.４㊀Ｃｏｎｖｅｒｇｅｎｃｅｏｆｎｅｕｒａｌｎｅｔｗｏｒｋｔｒａｉｎｉｎｇｅｒｒｏｒｓ采取同样的方式生成测试集来体现两种类型的神经网络代理模型计算效果ꎮ左侧流体的压力㊁密度在{０.２５ꎬ０.４５ꎬ０.６５}中取值ꎬ速度在{－２.０ꎬ－０.５ꎬ１.０ꎬ２.５}中取值ꎻ右侧流体的压力㊁密度在{０.３５ꎬ０.５５ꎬ０.７５}中取值ꎬ速度在{－２.５ꎬ－１.０ꎬ０.５ꎬ２.０}中取值ꎬ共计１２９６个测试样本ꎮ这些测试样本完全不同于训练样本ꎮ分别给出了两种类型的隐含层包含５个节点的神经网络模型和包含４０个节点的神经网络模型在该测试集上的误差值ꎬ如图５所示ꎮ横坐标表示样本点编号ꎬ纵坐标表示网络预测值与精确解的绝对误差ꎮ(ａ)Ｎｅｕｒａｌｎｅｔｗｏｒｋｗｉｔｈ５ｎｏｄｅｓｉｎｔｈｅｈｉｄｄｅｎｌａｙｅｒ(ｂ)Ｎｅｕｒａｌｎｅｔｗｏｒｋｗｉｔｈ４０ｎｏｄｅｓｉｎｔｈｅｈｉｄｄｅｎｌａｙｅｒ图５㊀网络预测值与精确解绝对误差对比图Ｆｉｇ.５㊀Ａｂｓｏｌｕｔｅｅｒｒｏｒｓｂｅｔｗｅｅｎｅｘａｃｔａｎｄｐｒｅｄｉｃｔｅｄｖｅｌｏｃｉｔｉｅｓ该结果初步表明ꎬ简化的神经网络模型的计算精度并没有因输入特征维度及样本量的减少而下降ꎬ并且每个隐含层包含４０个节点的神经网络模型的测试误差基本都在１％以下ꎮ３.２㊀一维气－气问题本节共选取两个一维空气－氦气问题测试各版本神经网络模型的计算效果ꎬ两个算例因初始状态的不同而产生不同的波系结构ꎮ８３第２期刘子岩ꎬ等:可压缩两气体流动的简化神经网络模型算例１的初始无量纲参数为ｘ<０.５:(γＬꎬｕＬꎬｐＬꎬρＬ)＝(１.４ꎬ－１.０ꎬ１.０ꎬ１.０)ｘ>０.５:(γＲꎬｕＲꎬｐＲꎬρＲ)＝(１.６６７ꎬ１.０ꎬ１.０ꎬ１.３８){在本算例中ꎬ物质界面两侧各产生一个稀疏波ꎬ运算到０.１５终止ꎮ算例２的初始无量纲参数为ｘ<０.４５:(γＬꎬｕＬꎬｐＬꎬρＬ)＝(１.４ꎬ－１.０ꎬ１.０ꎬ１.０)ｘ>０.４５:(γＲꎬｕＲꎬｐＲꎬρＲ)＝(１.６６７ꎬ１.０ꎬ１０.０ꎬ１.０){在本算例中ꎬ物质界面左侧会产生一个激波ꎬ界面右侧会产生一个稀疏波ꎬ运算到０.０８终止ꎮ分别采用基于隐式迭代的标准ＰＧＦＭ以及两种结构的神经网络代理模型的ＰＧＦＭ进行求解ꎬ在０<ｘ<１的计算域内共划分１０００个计算网格ꎮ空间离散使用５阶ＷＥＮＯ格式ꎬ时间离散使用３阶Ｒｕｎｇｅ￣Ｋｕｔｔａ格式ꎮ两个算例的速度㊁压力㊁密度计算结果对比分别见图６和图７ꎮ图中ＩｍｐｌｉｃｉｔＲＳ表示隐式迭代求解Ｒｉｅｍａｎｎ问题的方法ꎬ从图中可以看出ꎬ除５￣５￣５￣２㊁５￣１０￣１０￣２与３￣５￣５￣２等含较少神经元的神经网络出现了相对大的误差外ꎬ其他各网络均能较好地拟合标准ＰＧＦＭ结果ꎬ３￣４０￣４０￣２神经网络的计算结果几乎与标准ＰＧＦＭ计算结果完全重合ꎮ(ａ)Ｖｅｌｏｃｉｔｙ(ｂ)Ｐｒｅｓｓｕｒｅ(ｃ)Ｄｅｎｓｉｔｙ图６㊀算例１计算结果对比Ｆｉｇ.６㊀Ｃｏｍｐａｒｉｓｏｎｏｆｃａｌｃｕｌａｔｉｏｎｒｅｓｕｌｔｓｆｏｒｃａｓｅ１(ａ)Ｖｅｌｏｃｉｔｙ(ｂ)Ｐｒｅｓｓｕｒｅ(ｃ)Ｄｅｎｓｉｔｙ图７㊀算例２计算结果对比Ｆｉｇ.７㊀Ｃｏｍｐａｒｉｓｏｎｏｆｃａｌｃｕｌａｔｉｏｎｒｅｓｕｌｔｓｆｏｒｃａｓｅ２９３气体物理２０２４年㊀第９卷将算例２的界面处流体压力㊁速度和密度的部分神经网络预测结果与Ｒｉｅｍａｎｎ问题精确解进行对比ꎬ对相对误差定量统计ꎬ如表２所示ꎮ根据定量分析结果ꎬ简化的网络模型的误差甚至更小ꎮ表２㊀算例２相对误差定量统计表Ｔａｂｌｅ２㊀Ｑｕａｎｔｉｔａｔｉｖｅｒｅｓｕｌｔｓｏｆｒｅｌａｔｉｖｅｅｒｒｏｒｓｉｎｃａｓｅ２ｐｈｙｓｉｃａｌｑｕａｎｔｉｔｙｓｏｌｕｔｉｏｎ５￣１０￣１０￣２３￣１０￣１０￣２５￣４０￣４０￣２３￣４０￣４０￣２ｐＩｅｘａｃｔｓｏｌｕｔｉｏｎ２.６１２６５ｐｒｅｄｉｃｔｅｄｓｏｌｕｔｉｏｎ２.６９６１９２.６３０７９２.６０２５６２.６１３６１ｒｅｌａｔｉｖｅｅｒｒｏｒ３.１９８％０.６９４％０.３８６％０.０３７％ｕＩｅｘａｃｔｓｏｌｕｔｉｏｎ－１.８８３０４ｐｒｅｄｉｃｔｅｄｓｏｌｕｔｉｏｎ－１.８２２５４－１.８７０１５－１.８９０３５－１.８８２４２ｒｅｌａｔｉｖｅｅｒｒｏｒ３.２１３％０.６８５％０.３８８％０.０３３％ρＩＬｅｘａｃｔｓｏｌｕｔｉｏｎ１.９３６２１ｐｒｅｄｉｃｔｅｄｓｏｌｕｔｉｏｎ１.８５４０２１.９１２８４１.９４０８６１.９２９６８ｒｅｌａｔｉｖｅｅｒｒｏｒ４.２４５％１.２０７％０.２４０％０.３３７％ρＩＲｅｘａｃｔｓｏｌｕｔｉｏｎ０.４４７０１ｐｒｅｄｉｃｔｅｄｓｏｌｕｔｉｏｎ０.４５３５１０.４４６２６０.４４３８７０.４４５０１ｒｅｌａｔｉｖｅｅｒｒｏｒ１.４５４％０.１６８％０.７０２％０.４４７％可见在计算一维问题时ꎬ简化的神经网络模型也能够在训练效率和计算效率领先的情况下和５个输入特征的神经网络模型保持相近的计算精度ꎮ３.３㊀二维气－气问题本算例模拟左行激波(Ｍａｓ＝１.２２)扫过氦气气泡的过程ꎬ计算域如图８所示ꎬ该算例常被用于测试计算格式的性能ꎬ相关研究可参考文献[２０]ꎮ计算域上下边界采用反射边界条件ꎬ左右边界采用无反射边界条件ꎬ网格密度为６５１ˑ１７９ꎬＣＦＬ数设定为０.５ꎮ初始无量纲条件和状态方程相关参数为ａｉｒ(Ｉ):(γＩꎬｕＩꎬｖＩꎬｐＩꎬρＩ)＝(１.４ꎬ０.０ꎬ０.０ꎬ１.０ꎬ１.０)ａｉｒ(ＩＩ):(γＩＩꎬｕＩＩꎬｖＩＩꎬｐＩＩꎬρＩＩ)＝(１.４ꎬ－０.３９４ꎬ０.０ꎬ１.５６９８ꎬ１.３７６４)ｈｅｌｉｕｍｂｕｂｂｌｅ(ＩＩＩ):(γＩＩＩꎬｕＩＩＩꎬｖＩＩＩꎬｐＩＩＩꎬρＩＩＩ)＝(１.６６７ꎬ０.０ꎬ０.０ꎬ１.０ꎬ０.１３８)ìîíïïïïïïïï图８㊀空气中激波扫过氦气泡算例的计算域(单位:ｍｍ)Ｆｉｇ.８㊀Ｃｏｍｐｕｔａｔｉｏｎａｌｄｏｍａｉｎｏｆａｎａｉｒｓｈｏｃｋｉｍｐｉｎｇｉｎｇｏｎａｈｅｌｉｕｍｂｕｂｂｌｅ(ｕｎｉｔ:ｍｍ)使用嵌入物理约束的４０￣４０规模的神经网络模型对本算例进行计算ꎬ空间离散使用５阶ＷＥＮＯ格式ꎬ时间离散使用３阶Ｒｕｎｇｅ￣Ｋｕｔｔａ格式ꎮ同时ꎬ也给出基于隐式迭代的标准方法计算结果作为对比ꎮ图９给出３种方法在第７３１ꎬ１２０７ꎬ１６０４ꎬ３９９８时间步的流场结果ꎬ图９(ａ)为５￣４０￣４０￣２神经网络代理模型与基于隐式迭代的标准ＰＧＦＭ的对比纹影图ꎬ图９(ｂ)为３￣４０￣４０￣２神经网络代理模型与基于隐式迭代的标准ＰＧＦＭ的对比纹影图ꎮ从图中可以看出ꎬ激波从右向左扫过氦气泡ꎬ氦气泡在激波的作用下被压缩变形ꎬ并形成向四周扩散的膨胀波ꎮ膨胀波在上下壁面发生反射后进一步与气泡作用ꎬ最终形成复杂的波系结构ꎮ(ａ)５￣４０￣４０￣２ＮＮ＆ｉｍｐｌｉｃｉｔＲＳ０４第２期刘子岩ꎬ等:可压缩两气体流动的简化神经网络模型(ｂ)３￣４０￣４０￣２ＮＮ＆ｉｍｐｌｉｃｉｔＲＳ图９㊀激波扫过氦气泡算例纹影对比图Ｆｉｇ.９㊀Ｃｏｍｐａｒｉｓｏｎｏｆｔｈｅｓｃｈｌｉｅｒｅｎｆｏｒａｎａｉｒｓｈｏｃｋｉｍｐｉｎｇｉｎｇｏｎａｈｅｌｉｕｍｂｕｂｂｌｅ分析两个对比图ꎬ两种结构的神经网络与标准ＰＧＦＭ方法相比ꎬ均能清晰地反映波系结构ꎬ并不会产生明显的误差ꎮ可见在应用于二维问题时ꎬ简化的神经网络也能取得与５个输入量的神经网络相同的计算效果ꎮ４㊀结论本文基于之前提出的嵌入物理约束的可压缩多介质流神经网络模型ꎬ给出一种两气体流动的神经网络模型简化方法ꎮ首先详细介绍了在完全气体状态方程下输入特征从无界集到有界集的标准化映射方法ꎬ可提高神经网络模型的泛化能力ꎮ基于该标准化方法ꎬ通过选取较大压力㊁密度值作为参考压力㊁密度ꎬ可将神经网络的输入量从５个减少到３个ꎮ经过算例验证证实该方法的计算效果与５个输入量的神经网络模型相当ꎬ但神经网络的训练效果得到很大提高ꎮ另外ꎬ因为网络结构的简单化ꎬ计算效率也会进一步提高ꎬ该方法更加具备进一步开发的潜力ꎮ值得注意的是ꎬ目前该工作基于提到的完全气体状态方程具备的性质ꎬ因此目前只适用于两气体界面问题ꎮ对于含有其他状态介质的问题ꎬ还须继续展开相关研究ꎮ参考文献(Ｒｅｆｅｒｅｎｃｅｓ)[１]㊀刘铁钢ꎬ许亮.模拟多介质界面问题的虚拟流体方法综述[Ｊ].气体物理ꎬ２０１９ꎬ４(２):１￣１６.ＬｉｕＴＧꎬＸｕＬ.Ａｒｅｖｉｅｗｏｆｇｈｏｓｔｆｌｕｉｄｍｅｔｈｏｄｓｆｏｒｍｕｌｔｉ￣ｍｅｄｉｕｍｉｎｔｅｒｆａｃｅｓｉｍｕｌａｔｉｏｎ[Ｊ].ＰｈｙｓｉｃｓｏｆＧａｓｅｓꎬ２０１９ꎬ４(２):１￣１６(ｉｎＣｈｉｎｅｓｅ).[２]ＫａｗａｋａｍｉＫꎬＫｉｔａＹꎬＹａｍａｍｏｔｏＹ.Ｆｒｏｎｔ￣ｔｒａｃｋｉｎｇｓｉｍｕ￣ｌａｔｉｏｎｏｆｔｈｅｗｅｔｔｉｎｇｂｅｈａｖｉｏｒｏｆａｎｉｍｐｉｎｇｉｎｇｄｒｏｐｌｅｔｕｓｉｎｇａｒｅｌａｘｅｄｉｍｐｅｒｍｅａｂｉｌｉｔｙｃｏｎｄｉｔｉｏｎａｎｄａｇｅｎｅｒａｌｉｚｅｄＮａｖｉｅｒｂｏｕｎｄａｒｙｃｏｎｄｉｔｉｏｎ[Ｊ].Ｃｏｍｐｕｔｅｒｓ＆Ｆｌｕｉｄｓꎬ２０２３ꎬ２５１:１０５７３９.[３]ＯｓｈｅｒＳꎬＳｅｔｈｉａｎＪＡ.Ｆｒｏｎｔｓｐｒｏｐａｇａｔｉｎｇｗｉｔｈｃｕｒｖａｔｕｒｅ￣ｄｅｐｅｎｄｅｎｔｓｐｅｅｄ:ａｌｇｏｒｉｔｈｍｓｂａｓｅｄｏｎＨａｍｉｌｔｏｎ￣Ｊａｃｏｂｉｆｏｒｍｕｌａｔｉｏｎｓ[Ｊ].ＪｏｕｒｎａｌｏｆＣｏｍｐｕｔａｔｉｏｎａｌＰｈｙｓｉｃｓꎬ１９８８ꎬ７９(１):１２￣４９.[４]ＬｉｕＴＧꎬＫｈｏｏＢＣꎬＹｅｏＫＳ.Ｔｈｅｓｉｍｕｌａｔｉｏｎｏｆｃｏｍ￣ｐｒｅｓｓｉｂｌｅｍｕｌｔｉ￣ｍｅｄｉｕｍｆｌｏｗ.Ｉ.Ａｎｅｗｍｅｔｈｏｄｏｌｏｇｙｗｉｔｈｔｅｓｔａｐｐｌｉｃａｔｉｏｎｓｔｏ１Ｄｇａｓ￣ｇａｓａｎｄｇａｓ￣ｗａｔｅｒｃａｓｅｓ[Ｊ].Ｃｏｍｐｕｔｅｒｓ＆Ｆｌｕｉｄｓꎬ２００１ꎬ３０(３):２９１￣３１４.[５]ＬｉｕＴＧꎬＫｈｏｏＢＣꎬＹｅｏＫＳ.Ｇｈｏｓｔｆｌｕｉｄｍｅｔｈｏｄｆｏｒｓｔｒｏｎｇｓｈｏｃｋｉｍｐａｃｔｉｎｇｏｎｍａｔｅｒｉａｌｉｎｔｅｒｆａｃｅ[Ｊ].ＪｏｕｒｎａｌｏｆＣｏｍｐｕｔａｔｉｏｎａｌＰｈｙｓｉｃｓꎬ２００３ꎬ１９０(２):６５１￣６８１.[６]许亮ꎬ冯成亮ꎬ刘铁钢.虚拟流体方法的设计原则[Ｊ].计算物理ꎬ２０１６ꎬ３３(６):６７１￣６８０.ＸｕＬꎬＦｅｎｇＣＬꎬＬｉｕＴＧ.Ｄｅｓｉｇｎｐｒｉｎｃｉｐｌｅｓｏｆｇｈｏｓｔｆｌｕｉｄｍｅｔｈｏｄ[Ｊ].ＣｈｉｎｅｓｅＪｏｕｒｎａｌｏｆＣｏｍｐｕｔａｔｉｏｎａｌＰｈｙ￣ｓｉｃｓꎬ２０１６ꎬ３３(６):６７１￣６８０(ｉｎＣｈｉｎｅｓｅ).[７]ＸｕＬꎬＦｅｎｇＣＬꎬＬｉｕＴＧ.Ｐｒａｃｔｉｃａｌｔｅｃｈｎｉｑｕｅｓｉｎｇｈｏｓｔｆｌｕｉｄｍｅｔｈｏｄｆｏｒｃｏｍｐｒｅｓｓｉｂｌｅｍｕｌｔｉ￣ｍｅｄｉｕｍｆｌｏｗｓ[Ｊ].ＣｏｍｍｕｎｉｃａｔｉｏｎｓｉｎＣｏｍｐｕｔａｔｉｏｎａｌＰｈｙｓｉｃｓꎬ２０１６ꎬ２０(３):６１９￣６５９.[８]ＸｕＬꎬＹａｎｇＷＢꎬＬｉｕＴＧ.Ａｎｉｎｔｅｒｆａｃｅｔｒｅａｔｍｅｎｔｆｏｒｔｗｏ￣ｍａｔｅｒｉａｌｍｕｌｔｉ￣ｓｐｅｃｉｅｓｆｌｏｗｓｉｎｖｏｌｖｉｎｇｔｈｅｒｍａｌｌｙｐｅｒｆｅｃｔｇａｓｅｓｗｉｔｈｃｈｅｍｉｃａｌｒｅａｃｔｉｏｎｓ[Ｊ].ＪｏｕｒｎａｌｏｆＣｏｍｐｕｔａｔｉｏｎａｌＰｈｙｓｉｃｓꎬ２０２２ꎬ４４８:１１０７０７.[９]刘子岩ꎬ许亮.嵌入物理约束的可压缩多介质流神经网络模型[Ｊ].计算物理ꎬ２０２３ꎬ４０(６):７６１￣７６９.ＬｉｕＺＹꎬＸｕＬ.Ｎｅｕｒａｌｎｅｔｗｏｒｋｍｏｄｅｌｓｏｆｃｏｍｐｒｅｓｓｉｂｌｅｍｕｌｔｉ￣ｍｅｄｉｕｍｆｌｏｗｓｅｍｂｅｄｄｅｄｗｉｔｈｐｈｙｓｉｃａｌｃｏｎｓｔｒａｉｎｔｓ[Ｊ].ＣｈｉｎｅｓｅＪｏｕｒｎａｌｏｆＣｏｍｐｕｔａｔｉｏｎａｌＰｈｙｓｉｃｓꎬ２０２３ꎬ４０(６):７６１￣７６９(ｉｎＣｈｉｎｅｓｅ).[１０]ＸｕＬꎬＬｉｕＴＧ.Ｅｘｐｌｉｃｉｔｉｎｔｅｒｆａｃｅｔｒｅａｔｍｅｎｔｓｆｏｒｃｏｍｐｒｅｓ￣ｓｉｂｌｅｇａｓ￣ｌｉｑｕｉｄｓｉｍｕｌａｔｉｏｎｓ[Ｊ].Ｃｏｍｐｕｔｅｒｓ＆Ｆｌｕｉｄｓꎬ２０１７ꎬ１５３:３４￣４８.[１１]ＲａｙＤꎬＨｅｓｔｈａｖｅｎＪＳ.Ａｎａｒｔｉｆｉｃｉａｌｎｅｕｒａｌｎｅｔｗｏｒｋａｓａｔｒｏｕｂｌｅｄ￣ｃｅｌｌｉｎｄｉｃａｔｏｒ[Ｊ].ＪｏｕｒｎａｌｏｆＣｏｍｐｕｔａｔｉｏｎａｌＰｈｙ￣ｓｉｃｓꎬ２０１８ꎬ３６７:１６６￣１９１.[１２]ＦｅｎｇＹＷꎬＬｉｕＴＧꎬＷａｎｇＫ.Ａｃｈａｒａｃｔｅｒｉｓｔｉｃ￣ｆｅａｔｕｒｅｄｓｈｏｃｋｗａｖｅｉｎｄｉｃａｔｏｒｆｏｒｃｏｎｓｅｒｖａｔｉｏｎｌａｗｓｂａｓｅｄｏｎｔｒａｉｎ￣ｉｎｇａｎａｒｔｉｆｉｃｉａｌｎｅｕｒｏｎ[Ｊ].ＪｏｕｒｎａｌｏｆＳｃｉｅｎｔｉｆｉｃＣｏｍｐｕ￣ｔｉｎｇꎬ２０２０ꎬ８３(１):２１.１４气体物理２０２４年㊀第９卷[１３]ＦｅｎｇＹＷꎬＬｉｕＴＧ.Ａｃｈａｒａｃｔｅｒｉｓｔｉｃ￣ｆｅａｔｕｒｅｄｓｈｏｃｋｗａｖｅｉｎｄｉｃａｔｏｒｏｎｕｎｓｔｒｕｃｔｕｒｅｄｇｒｉｄｓｂａｓｅｄｏｎｔｒａｉｎｉｎｇａｎａｒｔｉｆｉ￣ｃｉａｌｎｅｕｒｏｎ[Ｊ].ＪｏｕｒｎａｌｏｆＣｏｍｐｕｔａｔｉｏｎａｌＰｈｙｓｉｃｓꎬ２０２１ꎬ４４３:１１０４４６.[１４]丘润荻ꎬ王静竹ꎬ黄仁芳ꎬ等.改进的物理融合神经网络在瑞利￣泰勒不稳定性问题中的应用[Ｊ].力学学报ꎬ２０２２ꎬ５４(８):２２２４￣２２３４.ＱｉｕＲＤꎬＷａｎｇＪＺꎬＨｕａｎｇＲＦꎬｅｔａｌ.Ｔｈｅａｐｐｌｉｃａｔｉｏｎｏｆｍｏｄｉｆｉｅｄｐｈｙｓｉｃｓ￣ｉｎｆｏｒｍｅｄｎｅｕｒａｌｎｅｔｗｏｒｋｓｉｎＲａｙｌｅｉｇｈ￣Ｔａｙｌｏｒｉｎｓｔａｂｉｌｉｔｙ[Ｊ].ＣｈｉｎｅｓｅＪｏｕｒｎａｌｏｆＴｈｅｏｒｅｔｉｃａｌａｎｄＡｐｐｌｉｅｄＭｅｃｈａｎｉｃｓꎬ２０２２ꎬ５４(８):２２２４￣２２３４(ｉｎＣｈｉ￣ｎｅｓｅ).[１５]ＭａｇｉｅｒａＪꎬＲａｙＤꎬＨｅｓｔｈａｖｅｎＪＳꎬｅｔａｌ.Ｃｏｎｓｔｒａｉｎｔ￣ａｗａｒｅｎｅｕｒａｌｎｅｔｗｏｒｋｓｆｏｒＲｉｅｍａｎｎｐｒｏｂｌｅｍｓ[Ｊ].ＪｏｕｒｎａｌｏｆＣｏｍｐｕｔａｔｉｏｎａｌＰｈｙｓｉｃｓꎬ２０２０ꎬ４０９:１０９３４５. [１６]ＷａｎｇＪＣꎬＨｉｃｋｅｙＪＰ.ＦｌｕｘＮｅｔ:ａｐｈｙｓｉｃｓ￣ｉｎｆｏｒｍｅｄｌｅａｒｎｉｎｇ￣ｂａｓｅｄＲｉｅｍａｎｎｓｏｌｖｅｒｆｏｒｔｒａｎｓｃｒｉｔｉｃａｌｆｌｏｗｓｗｉｔｈｎｏｎ￣ｉｄｅａｌｔｈｅｒｍｏｄｙｎａｍｉｃｓ[Ｊ].ＣｏｍｐｕｔｅｒＭｅｔｈｏｄｓｉｎＡｐ￣ｐｌｉｅｄＭｅｃｈａｎｉｃｓａｎｄＥｎｇｉｎｅｅｒｉｎｇꎬ２０２３ꎬ４１１:１１６０７０. [１７]ＲｕｇｇｅｒｉＭꎬＲｏｙＩꎬＭｕｅｔｅｒｔｈｉｅｓＭＪꎬｅｔａｌ.Ｎｅｕｒａｌ￣ｎｅｔ￣ｗｏｒｋ￣ｂａｓｅｄＲｉｅｍａｎｎｓｏｌｖｅｒｆｏｒｒｅａｌｆｌｕｉｄｓａｎｄｈｉｇｈｅｘｐｌｏ￣ｓｉｖｅｓꎻａｐｐｌｉｃａｔｉｏｎｔｏｃｏｍｐｕｔａｔｉｏｎａｌｆｌｕｉｄｄｙｎａｍｉｃｓ[Ｊ].ＰｈｙｓｉｃｓｏｆＦｌｕｉｄｓꎬ２０２２ꎬ３４(１１):１１６１２１.[１８]ＸｕＬꎬＬｉｕＺＹꎬＦｅｎｇＹＷꎬｅｔａｌ.Ｐｈｙｓｉｃｓ￣ｃｏｎｓｔｒａｉｎｅｄｎｅｕｒａｌｎｅｔｗｏｒｋｓａｓｍｕｌｔｉ￣ｍａｔｅｒｉａｌＲｉｅｍａｎｎｓｏｌｖｅｒｓｆｏｒｃｏｍｐｒｅｓｓｉｂｌｅｔｗｏ￣ｇａｓｓｉｍｕｌａｔｉｏｎｓ[Ｊ].ＪｏｕｒｎａｌｏｆＣｏｍｐｕ￣ｔａｔｉｏｎａｌＳｃｉｅｎｃｅꎬ２０２４ꎬ７８:１０２２６１.[１９]ＴｏｒｏＥＦ.Ｒｉｅｍａｎｎｓｏｌｖｅｒｓａｎｄｎｕｍｅｒｉｃａｌｍｅｔｈｏｄｓｆｏｒｆｌｕｉｄｄｙｎａｍｉｃｓ:ａｐｒａｃｔｉｃａｌｉｎｔｒｏｄｕｃｔｉｏｎ[Ｍ].３ｒｄｅｄ.Ｂｅｒ￣ｌｉｎꎬＨｅｉｄｅｌｂｅｒｇ:Ｓｐｒｉｎｇｅｒꎬ２００９:１５１￣１６２.[２０]ＣｏｒａｌｉｃＶꎬＣｏｌｏｎｉｕｓＴ.Ｆｉｎｉｔｅ￣ｖｏｌｕｍｅＷＥＮＯｓｃｈｅｍｅｆｏｒｖｉｓｃｏｕｓｃｏｍｐｒｅｓｓｉｂｌｅｍｕｌｔｉｃｏｍｐｏｎｅｎｔｆｌｏｗｓ[Ｊ].ＪｏｕｒｎａｌｏｆＣｏｍｐｕｔａｔｉｏｎａｌＰｈｙｓｉｃｓꎬ２０１４ꎬ２７４:９５￣１２１.２４。

honours 数列-回复数列（Sequences）是数学中一个非常重要且广泛应用的概念。

数列指的是按照特定规律排列的一组数。

在数学中，我们常常使用数列来表示各种各样的模式和趋势，从而帮助我们更好地理解和分析问题。

本文将以“数列”为主题，逐步回答与该主题相关的问题。

数列的定义是什么？数列是按照一定规律排列的一组数。

它通常用一对大括号或中括号来表示。

数列中的每个数被称为该数列的项（term）。

项的顺序通常用小标来表示，例如，第一个项用a1表示，第二个项用a2表示，依此类推。

数列的通项公式（general term formula）用来表示数列中第n项的表达式。

数列有哪些分类？数列可以被划分为不同的类型，其中最常见的几种包括等差数列、等比数列、斐波那契数列和调和数列。

等差数列（Arithmetic Sequence）是一种数列，每一项与它前面的项之间的差值都相等。

例如，1, 3, 5, 7, 9就是一个等差数列，其中每一项与前一项的差值都为2。

等比数列（Geometric Sequence）是一种数列，每一项与它前面的项之间的比值都相等。

例如，2, 4, 8, 16, 32就是一个等比数列，其中每一项与前一项的比值都为2。

斐波那契数列（Fibonacci Sequence）斐波那契数列是一个特殊的数列，从第三项开始，每一项都是前两项之和。

例如，1, 1, 2, 3, 5, 8, 13, 21就是一个斐波那契数列，其中第三项为前两项的和，第四项为前两项的和，依此类推。

调和数列（Harmonic Sequence）调和数列是一种特殊的数列，其中每一项都是其所有前一项的倒数之和。

例如，1, 1/2, 1/3, 1/4, 1/5就是一个调和数列，其中第二项为前一项的倒数，第三项为前两项倒数之和，以此类推。

数列有什么应用？数列是数学中非常重要的工具，它在很多领域都有广泛的应用。

首先，数列在数学分析和证明中扮演着重要的角色。

CompuTherm Training ClassShanghai, ChinaNovember 15, 2007CompuTherm, LLCMadison, WI 53179, USACopyright © 2007 CompuTherm LLCCompuTherm Training Class, 2007CompuTherm Training Class, 200711.1 OptimizationThis material is intended for training users on how to use PanOptimizer with Al-Zn as an example. Understanding optimization process itself is more complicated and requires extensive experience.Users should follow the seven steps below to perform the model parameter optimization using PanOptimizer. Here we use Al-Zn binary as an example.9 Step 1: Load Thermodynamic Database FileLoad thermodynamic database AlZn_0.TDB file into Pandat workspace. This database file contains the phases to be modeled. The database file definition is consistent with the current accepted format in the CALPHAD society.9 Step 2: Define Model ParametersPanOptimizer asks user to define the optimization parameters in the TDB file. Each model parameter is defined by the keyword “OPTIMIZATION ”. The format of defining a model parameter is:Optimization [name of model parameter][low bound][initial value][high bound] N !For example:// Keyword Name Low Bound Init. Value High Bound Optimization LIQ_AA 0 0; 60000 N ! Optimization LIQ_AAT -20 0; 20 N !Phase LIQUID % 1 1 !Constituent LIQUID : AL, ZN : !Parameter G(LIQUID, AL; 0) 298.15 G_Al_LIQUID; 6000 N ! Parameter G(LIQUID, ZN; 0) 298.15 G_Zn_LIQUID; 6000 N ! Parameter G(LIQUID, AL, ZN; 0) 298.15 LIQ_AA + LIQ_AAT*T; 6000 N !9 Step 3: Prepare Experimental FileThe most widely accepted format for experimental data file in the CALPHAD society is a POP file. PanOptimizer accepts most of the keywords in the POP format and adds a few special keywords. Typical format of a POP file to define a tie line looks like:CREATE_NEW_EQUILIBRIUM 1,1 CHANGE_STATUS PHASE * = SCHANGE_STATUS PHASE LIQUID, FCC_A1 = FIX 1 SET_CONDITION P = P0, T = 656EXPERIMENT X(LIQUID, AL)= 0.33 : DX1, EXPERIMENT X(FCC_A1, AL)= 0.11 : DX1,CompuTherm Training Class, 2007 29 Step 4: Load and Compile Experimental FileThrough the function Load Experimental File , AlZn_Sample.pop file can be loaded into the Pandat workspace and compiled with the database file loaded in step 1.9 Step 5: OptimizationModel parameter optimization is controlled through the optimization control panel as shown below. There are four control areas in the optimization control panel. They are: Histogram (A), Bound/Unbound Variables (B), Optimization (C), and Optimization Results (D). Histogram (A)Histogram is for displaying and tracing the history of the discrepancy between model-calculated values and experimental data. Free Bound/Unbound Variables (B)The user can choose the model parameters to be optimized in either the bounded or the unbounded mode. In the bounded mode, the low and highbounds defined in a TDB file will take into effects.CompuTherm Training Class, 20073Optimization (C)The optimization process can be controlled by choosing 1 Iteration or Run . In this example, there are totally 11 parameters to be optimized and all the initial values are set to be zero. The available experimental data include the enthalpy of liquid phase at 953K, two invariant reactions at 655K: LIQUID -> FCC_A1 + HCP_A3 and at 550K: FCC_A1#2 -> FCC_A1#1 + HCP_A3 and tie-lines: FCC_A1#1 + FCC_A1#2, LIQUID + FCC_A1, LIQUID + HCP_A3 and HCP_A3 + FCC_A1 at different temperatures.Before the optimization, user can check the calculated results using the initial model parameters without optimization. The parameters without optimization result in large discrepancies between the calculated values and the experimental measurements. Optimization of these model parameters (value zero) is needed.AB C DCompuTherm Training Class, 20074By clicking Run twice (two rounds), and each Run with 50 number of function calls, the sum of squares decreases from 904896 to 494.After another two rounds of Run, the sum of squares decreases from 494 to 2.8. Now, we can do the real time calculation using the instantly obtained optimized parameters. The following two comparisons show excellent agreements between the calculated results and the measurement data.CompuTherm Training Class, 2007Optimization Results (D)•During optimization, user can check model parameters through the Parameters button.•For each model parameter, user can change its low bound, upper bound and initial value in the dialog window.•It also allows user to Include or Exclude some model parameters for optimization using “Ctrl” or “Shift” to select those parameters.•If a set of optimized parameters is satisfactory, user can save this set of values as default ones through Set As Default. Otherwise, user can reject this set of values and go back to previous default values through Reset Default.5CompuTherm Training Class, 2007 6• User also can save the TDB file with the optimized model parameters through Save TDB .• To track the change of each parameter during the past 100 iterations, user can click the parameter name in the table.•It should be noted that any changes on model parameters made by user will take into effect only after the “Apply” button is clicked.Comparison with Experimental Data (D)• Through Experimental Data , users can compare the calculated values and the experimental data for each equilibrium.• The experimental information such as uncertainty and weight factor can be adjusted as needed during optimization procedure.•Through the column of “Residual ”, user can see different color after each comparison. The green means satisfactory results are achieved, while the red color means large deviation remaining.• User can Include or Exclude the selected experimental data for the optimization as needed.• More detailed phase equilibrium information can be accessed by clicking the equilibrium ID in the table.•The “Export Table” function is used to export the table data to Excel.CompuTherm Training Class, 200779 Step 6: Plot Diagram with Pptimized ParametersDuring optimization, user can calculate phase equilibria or other properties at anytime using the instantly obtained optimized parameters, and plot the diagrams with experimental data. This provides a graphic presentation on the comparison between the calculated and experimental results. The parameters used for calculation are displayed in the Parameters window.CompuTherm Training Class, 200789 Step 7: Save/Load Optimization ResultsDuring optimization, the user may save and load optimization file through Save Optimization Results and Load Optimization Results . The optimization results file has the extension name of POR that can only be read by PanOptimizer. All settings of the optimization are stored and recovered through these two operations. Examples on the settings are: what model parameters used in the optimization, and their changes from previous step; what experimental data used in the optimization, and their weight and uncertainty; and comparisons between calculated results and experimental data.1.2 Rough Search“Rough Search ” is used to obtain a rough topology of a phase diagram based on the given experimental data. It might be useful for quickly guessing the initial values of the model parameters at the beginning stage of theoptimization procedure.9Step 1: Load Thermodynamic Database File (same asSection 1.1)9Step 2: Define Model Parameters (same as Section 1.1)9Step 3: Prepare Experimental FileKeyword CREATE_ROUGH_EQUILIBRIUM is used to define phase equilibrium for rough search. When create a rough equilibrium, the status of phases in equilibrium can only be FIXED. The keyword PHASE_POINT is used to define the calculation conditions. No experimental data are allowed to be optimized for multiple-phase equilibria. However, a special case is the single-phase equilibrium. The thermodynamic property data of the single-phase equilibrium can be used as experimental data to optimize the model parameters. Example is:$ Tieline between HCP_A3 and FCC_A1 at 500K$ 10 is the weight factor of this equilibriumCREATE_ROUGH_EQUILIBRIUM @@, 1, 10CHANGE_STATUS PHASE FCC_A1, HCP_A3 = FIX 1PHASE_POINT P = P0, T = 500PHASE_POINT X(FCC_A1, AL)= 0.92, X(HCP_A3, AL)= 0.01For a phase described by a sublattice model, the initial value of an element or species in a sublattice should be specified using SET_START_VALUE.9Step 4: Load and Compile Experimental FileThrough the function Load Experimental File, AlZn_Sample_rough.pop file can be loaded into the Pandat workspace and compiled with the database file loaded in step 1. The “Rough Search” button in toolbar will be activated automatically if there are rough equilibria loaded into PanOptimizer.9Step 5: Rough SearchAgain, there are four control areas in the control panel for rough search as shown below (only Area C is marked since the other three areas are the same as the normal optimization control panel).CRough Search (C)There are two types of built-in optimization algorithms for Rough Search and they are Global Search and Local Search, as shown in area C.•Global Search is to search a set of good initial values of the model parameters within the whole searching domain, the initial values defined in TDB file will be ignored.•Local Search leads to a quick convergence to an optimal solution based on the given initial values of model parameters.The user can set the maximum number of function calls for both algorithms.For example, 1000 function evaluations of Global search lead to the sum of squares of 2.91. The figure below shows the calculated phase boundaries after a global search with 1000 function evaluations. It can be seen that the optimized parameters give a very similar topology with the known Al-Zn phase diagram. By clicking Local Search , the sum of squares decrease from 2.91 to 0.059, which reproduces the phase diagram in better agreement with the experimental one as you can see below.Optimization Results (D)•Model parameters can be checked through Parameters .•Chemical potential differences between any two phases in equilibrium can be accessed by Experimental Data . They are reflected by the values in column “(Wtd. Residual)^2” as shown in the dialog window below.• The weight factor of the equilibrium can be adjusted as needed during optimization procedure.• Chemical potential values of each phase in equilibrium can be accessed by clicking the buttons in the first column.• The “Export Table ” function is used to export the table data to Excel.9Step 6: Go Back to Normal OptimizationOnce the phase diagram topology is obtained, user can append the AlZn_sample.pop, and go through Steps 5-7 to obtain an optimal set of model parameters which can best fit phase diagram and thermodynamic properties.。