Julia和R软件在多介质动力学模型建模中的比较

莫俊超，吴孝槐，舒耀皋

（1. 上海化工研究院有限公司 检测中心，上海 200062；2. 上海化学品公共安全工程技术研究中心，上海 200062；3. 中国合格评定国家认可委员会，北京 100062）

当前，我国环境管理政策逐渐从末端治理转向风险管理［1-2］。多介质模型是评估化学物质环境风险的常用方法，按照模型的结构可分为一级、二级、三级和四级模型［3-4］。在实际环境系统中，化学物质的含量受人类活动排放、降解和迁移等过程的影响，常随时间发生改变，因此，能描述化学物质含量变化的多介质动力学模型（即四级模型）正得到越来越多的研究和应用。

多介质动力学模型的求解方法可分为数值解法［5］和解析解法。当多介质动力学模型中微分方程数量超过3个，或者模型中化学物质的排放速率、相体积、相流速和温度不恒定时，模型的解析解很难求出［6］，这时求解模型的数值解便成为了唯一的求解方法。

求解多介质动力学模型的数值解需要使用科学计算软件，目前使用较多的软件包括Matlab［7-10］、Mathematica［11］、C++［12］、Anylogic［13］、VB［14-15］、Python［16］和R［17］等。其中，Matlab、Mathematica、Anylogic和VB，以及部分C++编译器属于商业软件，价格昂贵，更新较慢，且对版权要求较多，限制了其应用范围。近年来，免费、开源且更新迅速的科学计算软件发展很快，其丰富的包和库等功能让用户可使用前人已有的成果，不必从零开始编写计算程序，从而大大减少了编程的工作量。该类软件以Python和R为代表，但目前未见系统讨论其在多介质动力学模型建模方面应用的研究报道。

Julia是近10年内开发的一种主要面向科学计算的软件，首个版本发布于2013年。同Python和R一样，Julia也是一款自由软件，并且结合了已有科学计算软件的优点。目前未见应用Julia进行多介质动力学模型建模的报道。

考虑到Python中算法的劣势［18］，本研究选择Julia和R，通过构建3个不同复杂程度的多介质动力学模型，对Julia和R在程序编写、运行时间、算法选择、敏感性分析和便捷性等方面进行比较，得出针对不同的多介质动力学模型适用的软件，以期拓展Julia和R在多介质动力学模型建模中的应用。

1 材料与方法

1.1 软件和包

使用R 4.0.2和Julia 1.5.2进行模型构建。R的下载地址为https：//cran.r-project.org/，使用的包有deSolve（版本1.28）、ODEsensitivity（版本1.1.2）和diffeqr（版本1.0.0）；Julia的下载地址为https：//julialang.org/，使用的包有DifferentialEquations.jl［19］（版本6.15.0）、LSODA.jl（版本0.6.2）、Sundials.jl（版本4.3.0）和DiffEqSensitivity.jl（版本6.31.6）。

1.2 EQC模型简介

选择一个通用模型EQC（equilibrium criterion，平衡判据）和EQC模型中的标准环境进行讨论。EQC模型于1996年被首次提出［20］，现已广泛应用于化学物质的风险筛查和评估中［21-25］。EQC模型中包括水、大气、土壤和沉积物4个环境介质，EQC动力学模型即对每个环境介质中化学物质的物质的量编写微分平衡方程：dM/dt=v输入-v输出，式中M为化学物质的物质的量，t为时间，v输入和v输出分别为输入速率和输出速率。其构建过程详见参考文献［26］。

1.3 不同复杂程度动力学模型的构建

为了构建不同复杂程度的动力学模型，选择不同函数作为EQC动力学模型中的排放速率函数，同时将模型的结束时间延长，具体参数见表1。一般来说，化学物质不会直接排放进入沉积物中，因此将沉积物中的排放速率设为0。所有模型中步长均设为0.1 h。

表1 不同复杂程度的模型参数

1.4 物质选择

为了对比不同性质的物质在模型中的运行情况，选择4种性质相差较大的物质进行模拟，其参数详见表2［26-27］，表中PCB-209为十氯联苯，KOW为正辛醇-水分配系数。

1.5 函数、算法和计算机配置

R中，使用deSolve包中的ode函数进行模型计算，计算方法选择ode函数中的lsoda（即ode函数的默认算法）、bdf和radau，使用system.time函数获取程序的运行时间。Julia中，使用DifferentialEquations.jl包中的solve函数进行模型计算，计算方法选择lsoda（）（即对应R中ode函数默认算法lsoda的Julia算法）、AutoTsit5（Rosenbrock23（））和CVODE_BDF（），其中lsoda（）算法需要加载LSODA.jl包，CVODE_BDF（）算法需要加载Sundials.jl包。Julia中使用@time命令获取程序的运行时间。计算所用的计算机处理器为Intel®Core i7-8565U，内存为8 GB，操作系统为Microsoft®Windows 10专业版。

2 结果与讨论

2.1 程序编写

由于所编写的程序较复杂，不再列出。R的主体程序（即不包括参数赋值计算的语句）为14行，Julia的主体程序为16行，R较Julia程序更简洁。R程序只经过修改一些基本语法便可在Julia上运行，反之亦然，两者语法非常类似，可移植性较强。

2.2 模型运行时间

对于每个程序，运行5次，计算得出运行时间的平均值，详见表3。

表3 R和Julia中EQC动力学模型的运行时间

2.3 模型敏感性分析方法和运行时间

模型的敏感性分析方法主要分为全局敏感性分析和局部敏感性分析两个方面，具体可参考陈卫平等［28］的总结。在R和Julia中，全局敏感性分析方法均有已发布的包可使用，用户不必自行编写敏感性分析程序，使用非常方便。R中模型全局敏感性分析包ODEsensitivity，提供了Morris法和Sobol法两种敏感性分析方法；Julia中DiffEqSensitivity.jl包可用来进行全局敏感性分析，包含了Morris、Sobol、eFAST和Regression 4种方法。另外，Julia中还提供了局部敏感性分析的函数，而R中暂未见相关功能。使用R和Julia中Morris和Sobol方法对表2中4个半衰期的敏感性进行分析，每个参数的值浮动±10%并设为均匀分布，计算模型敏感性分析的运行时间。求解算法均选择lsoda，Sobol法中模拟次数选择1 000次，结果见表4。

表2 所选4种物质的性质参数

表4 R和Julia中EQC动力学模型敏感性分析的运行时间

2.4 软件的比较

2.4.1 运行时间

由表3可见，对于相同的模型，Julia的运行时间至多为R的约1/10，尤其在复杂模型的lsoda算法中运行时间约为R的1/100，计算效率非常高。R的主要应用领域为数据统计，在科学计算方面的计算效率并没有进行很好的优化；Julia主要的应用领域为科学计算，从本研究的结果可见，其在多介质动力学模型计算中的计算效率是R的至少10倍。在敏感性分析方面，由表4可见，敏感性分析比模型求解需要更长的运行时间，Julia的运行效率同样是R的至少10倍。另外，不同物质在同一模型中的运行时间差别很小，这表明，化学物质性质不影响模型的运行时间。

2.4.2 算法

目前R中ode函数提供了包括lsoda、bdf和radau等算法在内的17种算法，可计算刚性和非刚性模型，默认算法lsoda可自动选择刚性或非刚性求解器。多介质动力学模型中不同介质间化学物质含量往往相差很大，大多数情况下模型为刚性模型，这时不能使用rk、rk4或euler等求解非刚性模型的算法进行计算，故R中适合多介质动力学模型求解的算法不多。本研究中选择了效率最高的3种算法，对比来看，lsoda和bdf算法的运行时间最短，计算效率最高，是R中求解多介质动力学模型的首选算法，但两者求解复杂模型的运行时间均很长。

Julia中求解微分方程模型的算法非常多，可用来进行多介质动力学模型求解的算法超过10种，为用户提供了丰富的选择。除了本研究中使用的lsoda（）、AutoTsit5（Rosenbrock23（））和CVODE_BDF（） 3种算法外，还可以尝试QNDF（）、TRBDF2（）和RadauIIA（）等多种算法。由表3可见，lsoda（）算法的运行时间最少，均不超过5 s，是Julia中是求解多介质动力学模型的首选算法。在敏感性分析方法中，Julia同样提供了丰富的选择，而R中目前仅提供了两种方法。

2.4.3 便捷性

由以上讨论可知，从运行时间和算法方面比较，Julia相对R有很大优势。然而，Julia发展时间较短，很多软件暂未提供Julia的接口，其应用生态暂没有R完善；Julia的中文资料仍然较少，中文社区也相对较小；目前国内Julia服务器很少，连接不稳定，Julia的安装相对R更繁琐。综上，Julia的应用暂不如R便捷。

为了结合R的便捷性和Julia的计算优势，RACKAUCKAS［29］在2020年8月首次发布了R中的diffeqr包，可令用户在R中直接调用Julia的微分方程计算函数，R中的JuliaCall包也可直接调用Julia中的其他函数。但目前diffeqr和JuliaCall包暂不完善，在针对复杂模型进行快速高效的模型计算和分析时，Julia仍然是较好的选择。

2.4.4 与其他软件的比较

结合以上结果和参考文献［18］，对多介质动力学模型求解中常用的软件R、Julia、Matlab、SciPy、Mathematica、C++、Anylogic和VB，从是否免费、算法选择、计算效率、模型分析工具、便捷性等方面进行比较，详见表5。综合来看，Julia在多介质动力学模型求解中优势较大，对于复杂的模型，建议首先选用Julia；但由于Julia的便捷性不如R，对于简单的多介质动力学模型，R也是可以满足需求的，这时使用R更方便。

表5 多介质动力学模型求解中的常用软件比较

3 结论

a）分别使用R和Julia，编写了4种化学物质在EQC动力学模型中的求解程序和半衰期的敏感性分析程序，得到模型求解程序和敏感性分析程序的运行时间。结果表明，Julia的计算效率是R的至少10倍，在一些情况下Julia的计算效率甚至是R的近100倍。化学物质性质不影响模型的运行时间。lsoda是Julia中求解多介质动力学模型的首选算法，而lsoda和bdf是R中的首选算法。

b）对多介质动力学模型求解中常用的软件，从是否免费、算法选择、计算效率、模型分析工具、便捷性等方面进行了比较。对于简单模型，建议使用R来建模、求解和分析；对于复杂模型，则建议使用Julia。

c）目前我国对多介质动力学模型的研发投入较少，缺乏适用于我国本土环境的多介质动力学模型，结合本研究的结果，建议考虑使用开源软件Julia和R来进行我国多介质动力学模型的研发和应用工作。