技术干货 | MindSpore AI科学计算系列（三）：SciML分析

本篇是MindSpore AI科学计算系列(一)续篇。我们将在上面介绍PINNs和英伟达SimNet工具包的基础上，重点分析AI科学计算领域的另一个软件包SciML。

点击↓回顾链接跳转查看文章

技术干货 | MindSpore AI科学计算(1)：PINNS和Simnet分析

SciML软件包

SciML是由NumFOCUS组织赞助的AI科学计算开源软件包。该组织旨在开发和维护一套功能丰富、性能高的科学机器学习工具，拥有非常活跃的社区，聚集了大量优秀的开发者。SimNet不同的是，SciML主要是基于Julia语言开发。这种语言是根据高性能数值分析和科学计算的需要而设计的。SciML具有优异的并行和分布式计算能力，它可以方便地利用社区中许多成熟的高性能基本代码。

图一 SciML研究主题：领域感知、可解释性和鲁棒性

图二 SciML应用主题：科学推理和数据挖掘、机器学习增强建模、智能自动化和决策支持

SciML的远景目标就是生成一个统一科学计算与机器学习的工具包，它包含科学计算与机器学习相融合领域的全部内容(见图一、二)。在传统科学计算方法集成的基础上，工具包进一步集成AI该领域的前沿算法提供了高度模块化、高扩展性和高性能的生态系统，以促进科学计算领域的基础和应用问题的研究。软件包除了微分方程求解器等科学计算模块外，还为正向求解、反问题、自动模型发现和决策推理提供了一套方法论。

SciML主要模块

与SimNet/DeepXDE不同，SciML软件包并不以TensorFlow或者Pytorch等比较成熟AI训练框架是后端。充分利用它。 Julia语言在科学计算领域的优势，发展出的一整套科学计算与机器学习相融合的生态系统。该工具包含大量应用库的应用库，下面只列出核心组件的一部分。

图三 SciML主模块及功能介绍

SciML优势

SciML软件包是科学机器学习和微分方程建模工具开发套件的集合，提供了丰富的工具库，形成了连贯的生态系统。与其他工具包相比，SciML优点如下：

1.多功能、多平台支持

从下表的对比中可以看出，SciML生态系统是唯一微分方程的生态系统（Stiff ODEs）、微分代数方程（DAEs）、随机微分方程（SDEs）、延迟微分方程（DDEs）、稳定邻域、分布式和多线程计算，集成深度学习的微分方程求解器库 。这些模板求解器支持数组维度的数组、数据类型和各种浮点精度的数值操作，并支持基于GPU/MPI/多线程等计算形式及其他功能。

表一 比较不同工具包的特性支持

2、高性能

SciML是基于Julia在科学计算领域，语言开发具有高性能的自然优势(如下图4所示，不同工具包的性能比较)。此外，性能也是如此SciML社区中所有关于性能的问题都将被视为优先事项bug支持科学机器学习的高性能扩展。最后，SciML跟上计算硬件的发展，确保与最新的高性能计算工具兼容。SciML随着新硬件技术的发布，编译器后端将不断更新，以维护可以在所有主要平台上工作的工具链。

图四 解决刚性的不同工具Hires性能比较方程

相比于SUNDIALS CVODE(C )以及Fortran方法（Radau）, SciML微分方程求解器库具有明显的性能优势。

算法层面，SciML也致力于提高性能。科学模型通常包含局部结构，特别是在离散方程中，通常只有相邻网格点之间的物理耦合。这种局部性导致了程序结构中张表达的稀疏性。SciML一套相互关联的工具可以直接在模拟模型上生成数值求解器代码，从而实现显著的性能提升。

3.高可靠性和鲁棒性

当前，SciML所有软件包都经过了科学模拟和机器学习工具的常规和鲁棒测试。在基准测试方面，SciML设置自动更新以获得更好的性能回归跟踪，并不断增加MPI和GPU新的基准SciML Benchmarks（https://github.com/SciML/SciMLBenchmarks.jl）。

支持DifferentialEquations.jl微分方程求解器和Flux.jl深度学习库之间的连接DiffEqFlux.jl工具仍在不断开发中。将代理模型Surrogates.jl与DifferentialEquations.jl和Zygote.jl常规测试等待机器学习辅助工具，确保代理模型工具能够在各种形式的微分方程上进行训练，然后在深度学习堆栈中使用。正是这种互联性将使下一代SciML从而影响方法论的生产化 大科学和工业应用。

4.丰富的加速算法开发和研究工具库

SciML致力于促进数值模拟方法和科学机器学习的尖端研究。基于此目标开发的许多工具可以自动完成基准测试和新方法的验证过程，以确保其安全性和可靠性，并加快这些方法向出版物和用户的转换。SciML通过DiffEqDevTools.jl和RoototedTrees.jl该模块使社区能够更多地利用和开发其提供的大型包装算法套件来快速测试和部署新算法。模块等外部工具的包装器FEniCS.jl和SciPyDiffEq.jl，跨平台比较很容易。另外，SciML它还提供了科学机器学习中最新算法的先进实现，如使用神经网络自动解决符号描述中的微分方程，并生成无网格离散器。

案例介绍

SciML广泛应用于微分方程、反问题、参数估计、决策推理、通用训练和神经常微分方程，建立了许多科学机器学习软件生态系统的基准案例。本节只提取了两个简单的案例介绍SciML使用工具集的方法。

神经常微分方程

本小节以Lotka-Volterra例如，如何使用方程SciML正向解决方程以及如何将其嵌入神经网络。Lotka-Volterra该方程是一种经典的动态模型，用于描述生物系统中捕食者和猎物种群的增长。该方程作为二元一级非线性微分方程组，具有以下形式：

SciML的DifferentialEquations.jl该模块提供了自定义精度、残差、求解器算法等多功能选项。给定初始条件后，参考以下代码调用solve计算接口可以直接给出如下图5所示(a)所示的解。

using DifferentialEquations, Plots ## Setup ODE function lotka_volterra(du, u, p, t) x, y = u alpha, beta, delta, gama = p du[1] = dx = alpha * x – beta * x * y du[2] = dy = -delta * y   gama * x * y end u0 = [1.0, 1.0] tspan = (0.0, 10.0) p = [1.5, 1.0, 3.0, 1.0] prob = ODEProblem(lotka_volterra, u0, tspan, p) ## Verify ODE solution Sol = solve(prob, Tsit5()) plot(Sol)

图五 （a）SciML正向求解Lotka-Volterra方程计算结果； (b) 神经常微分方程网络输出结果

与残余网络等经典结构相比，神经微分方程是以连续的方式对数据进行建模，从而连续化神经网络层次。下面我们将上述常微分方程嵌入神经网络的框架中，然后继续使用标准的优化技术，比如ADAM优化它们的权重。如下代码所示，微分方程层diffeq_rd接收积分的参数p，把它放在上面prob在定义的微分方程中，用选定的参数(求解器、残差等)求解。为了确定最小化损失函数的参数p，Flux训练这样一个网络：即神经网络的前向只包括解常微分方程。当猎物数量远离1时，损失函数会增加惩罚，因此训练结果中发现的捕食者和猎物族数量接近常数1（如图所示(b)所示)。在这种情况下，只有一层求解ODEs事实上，网络结构可以添加到神经网络的任何位置。

## Build a neural network that sets the cost as the difference from the generated data and 1 p = [2.2, 1.0, 2.0, 0.4] p=Flux.Params([p]) # Initial parameter vector function predict_rd() # 1-layer neural network diffeq_rd(p, prob, Tsit5(), saveat=0.1)[1, :] end loss_rd() = sum(abs2, x-1 for x in predict_td()) # loss function ## Optimize the parameter so the ODE’s solution stays near 1 data = Iterators.repeated((), 100) opt = ADAM(0.1) Flux.train!(loss_rd, [p], data, opt) plot(solve(remake(prob,p=Flux.data(p)),Tsit5(),saveat=0.1),ylim=(0,6))

2、PINNs松方程strog>

using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux                     
@parameters x, y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
## 2D PDE
Eq = Dxx(u(x, y) + Dyy(u(x, y)) ~ -sin(pi * x) * sin(pi * y)
## Boundary conditions
bcs = [u(0, y) ~ 0.f0, u(1, y) ~ -sin(pi*1) * sin(pi*y), u(x, 0) ~ 0.f0, u(x, 1) ~ -sin(pi*x) * sin(pi*1)]
## Space and time domains
domains = [x ∈ IntervalDomain(0.0, 1.0), y ∈ IntervalDomain(0.0, 1.0)]
## Discretization
dx = 0.1
## Neural network
dim = 2
chain = FastChain(FastDense(dim, 16, Flux.σ), FastDense(16, 16, Flux.σ), FastDense(16, 1))
discret = PhysicsInformedNN(chain, GridTraining(dx)) # PINNs model
pde_system = PDESystem(eq, bcs, domains, [x, y], [u])
prob = discretize(pde_system, discret)
res = GalacticOptim.solve(prob, Optim.BFGS(); maxiters=1000)
phi = discret.phi