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來進行我國多介質動力學模型的研發和應用工作。