动态系统建模被各领域广泛应用，例如电动汽车，能源系统，航空航天。我们本文提到动态系统主要是被控对象，对被控对象进行建模是因为我们希望了解这个系统（被控对象）的物理特性以及接受一些外部输入（力，扭矩，电流等等）时会有什么样的动态响应，基于此从而可以更好的给出控制输入得到我们期望的系统的输出，以及理解系统的退化或最大化提升系统效率。
这些动态系统的行为是由多物理场复杂的交互作用决定的，因此系统行为和系统响应建模通常需要复杂的第一原理支撑，仿真时也需要大量的计算(例如有限元模型)。
这也是本文的出发点，提供数据驱动（主要介绍深度学习和系统辨识）的模型降阶（Reduced Order Modeling）提速的方法，通过数据得到具有一定保真度的数据模型，在捕捉到系统动态特性的同时也提升仿真速度。
本文中将涉及多个 demo，数据以及脚本文件，若您感兴趣进一步获取这些链接，可以在文末填写反馈问卷，获取这些链接。
动态系统
动态系统包含状态空间 S, 时间集 T 和一个映射（规则）来描述状态随时间的演变规则 R:S×T→S。例如给定一个时刻 t 的状态 st ,通过这个规则可以计算后面一个或几个时刻状态 st+1=R(st)，st+2=R(R(st)) 等等。动态系统通常可以用随时间变化的方程或方程组来描述。尤其对于连续时间系统，可以通过微分方程来表示。
我们先看一个简单的常微分系统（ODE）, [链接1]
其中 y(t) 是系统状态。例如：一个简单二自由度线性系统，
其中 A 是一个 2x2 的矩阵。初始条件 x0 = [2; 0], 可以通过求解 ODE 方程得到的相应的二维输出 x(t)，包含两个状态，时序和动态图如下：
图表 1 系统输出x(t):
x0 = [2; 0];
A = [-0.1 -1; 1 -0.1];
trueModel = @(t,y) A*y; % 定义系统函数，此处就是一个状态空间方程
numTimeSteps = 2000;
T = 15;
odeOptions = odeset(RelTol=1.e-7);
t = linspace(0, T, numTimeSteps);
[~, xTrain] = ode45(trueModel, t, x0, odeOptions);
后面（在介绍 Neural ODE 部分）我们会尝试利用这个系统的数据 xTrain 进行深度学习模型的训练来得到这个系统的数据代理模型（Surrogate Model），这种思路也可以同样用于复杂系统。
既然数据驱动，有很多机器学习和深度学习算法可以用，那动态系统建模有什么特殊性呢，不是都适用吗？
目前工程中已经用到很多稳态（静态）模型。例如在发动机排放标定，通过 DoE 试验时我们会将发动机维持在不同的稳态工况（转速恒定，扭矩恒定等等），通过试验数据建模得到用于标定的稳态数据模型[链接2]。
稳态工况下，对于方程（1）这样一个简单系统，其中 y(t)' 可以看作 0，到达平衡点，于是 y(t) 和 u(t) 关系恒定，不再在时间维度上与历史状态 y(t-1)，y(t-2) 等等有关，因此稳态模型针对稳态工况是非常准的。
而在瞬态工况下通常 y(t)' 非零，因此方程（1）在求解系统输出 y(t) 时不仅由当前时刻的输入 u(t) 决定，还取决于 y(t-1)，y(t-2)，u(t-1), u(t-2) 等等，这就是动态系统的特殊性，当前输出不仅依赖于当前的输入，还依赖于系统过去的行为（历史输入和历史输出）。我们在下一篇关于系统辨识的文章中会基于示例详细说明这一点。
不考虑动态系统，单纯从系统建模来说，通常有如下的两个方向：基于第一原理的和基于数据驱动的。
第一原理建模是领域工程师都比较熟悉的，例如可以使用 M 脚本语言，Simulink 或 Simscape 建模语言从物理原理进行系统模型的创建。
那什么情况下搭建系统会考虑使用或部分使用数据代理模型（Surrogate Model）？
物理系统原理比较复杂或者不够清晰，无法构建第一原理模型
数据获取相对简单
第一原理物理模型过于复杂，求解耗时，在控制开发时有时效要求，需要加速仿真计算
可以看到右半部分主要是基于数据驱动的建模手段，其中针对动态系统的建模主要是系统辨识和神经网络/深度学习。也是本系列两个方向。
本文介绍神经网络的几种用于动态系统建模的模型，下一篇文章会介绍系统辨识的几种模型。
前馈神经网络（Feedforward Neural Network）
前神经网络大家都不陌生，各神经元分层排列。如图2
图2 简单前馈神经网络示意图
每层神经元只接受上层输出，结果只传给下一层，没有反馈。稍微复杂点的如图3。
图3 squeeze net MATLAB 示例
前馈神经网络是相对于循环神经网络（Recurrent Neural Network）而言，后者具有反馈，后面我们也会介绍。
那么问题来了，前馈神经网络因为没有反馈，如何表达动态系统（Dynamic System）的时间状态依赖呢？
通常主要手段使用不同尺度的滑动窗口来构建衍生特征从而表征系统在时间上的动态。
电池 SoC 预测示例
我们以电池管理系统为例，通过使用深度学习来估计电池的荷电状态 SoC（State of Charge）。
本示例主要介绍 SoC 数据驱动的建模方法，在即使不清楚电池电化学模型以及物理非线性特性的情况下，依然可以进行 SoC 估计。我们可以通过实验室中收集到的实测数据进行一个前馈神经网络代理模型训练实现 SoC 估计[1]。
脚本和数据[链接3]
图表4 数据集预览和模型的5个输入1个输出
数据准备
我们利用实测数据通过预处理得到训练数据集 cdsTrain，我们看其中一条数据记录：
preview(cdsTrain)
ans = 1×2 cell array
{5×669956 double} {1×669956 double}
数据包含五个特征，分别是电压 V、电流 I、温度 T、滑动平均电压 V_avg、滑动平均电流 I _avg,其中后面两个衍生特征用于表征动态信息。
搭建模型
接下来我们构建神经网络模型，结构比较简单，三层全连接网络，一个输出的 Feedforward Neural Network模型。
layers =[sequenceInputLayer(numFeatures,"Normalization","zerocenter")
fullyConnectedLayer(numHiddenUnits)
tanhLayer % HyperbolicTangent
fullyConnectedLayer(numHiddenUnits)
leakyReluLayer(0.3) % 激活函数
fullyConnectedLayer(numResponses)
clippedReluLayer(1) % 激活函数
regressionLayer];
设置训练选项并训练
options =trainingOptions('adam', ... % Adam optimizer
'MaxEpochs',Epochs,...
'ExecutionEnvironment','cpu',...%可以选择GPU
'InitialLearnRate',InitialLearnRate, ...
'LearnRateSchedule','piecewise', ...
'LearnRateDropPeriod',LearnRateDropPeriod,...
'LearnRateDropFactor',LearnRateDropFactor,...
'ValidationData', {X,Y}, ...
'ValidationFrequency',validationFrequency,...
'MiniBatchSize',miniBatchSize, ...
'CheckpointPath', NET_Path);
图表 5 训练过程 Loss 变化
导入测试数据验证模型在测试集上的准确度
图表 6 深度学习模型的预测值与实测值比较
仿真测试以及代码生成
Simulink 中的深度学习推断模块[链接8]支持将我们训练好的模型作为 block 参数, 一起作为被控对象集成到整个电池管理系统中。
图表 7 电池管理系统和电池的系统模型
图中所示的 BMS 的 Simulink 模型可以监控电池状态，确保运行安全，还有一个电池模型用于仿真电池的动态和负荷。上面训练的深度学习 SoC 预测器和其他电池平衡逻辑 Block 一样嵌入在 BMS 中可以闭环仿真，以及后面做代码生成与硬件在环。
图表 8 Simulink 中原生的用于深度学习推断的 Block 和 BMS 中的闭环测试 SoC 预测效果
图表 9 deep learning 模块的 C 代码生成
上面这个示例我们利用衍生特征（通过时间滑窗构建时域依赖的特征，也可以通过不同长度的滑窗构建多尺度的时域特征），将前馈神经网络用于动态系统的建模。
Temporal convolutional network (TCN)
TCN 的主要构成是一个扩展因果卷积层。任何一个时刻的计算输出都是基于多个历史时刻的输入。
它用于构建动态系统的逻辑和前面介绍的衍生特征是类似的，都是考虑了历史的多个时间步的输入，因此也可以用于建模动态系统。
TCN 从之前的时间步构建依赖，通常需要将多个卷积层叠加在一起。为了获得更长期的依赖关系，卷积层的膨胀因子呈指数级增加，如下图所示。
假设第 k 个卷积层的膨胀因子为2⁽ᵏ⁻¹⁾ ,步长为 1，则该网络的考虑到的依赖的时间窗的大小可计算为 R=(f-1)(2ᵏ-1)+1,其中 f 为过滤器大小，K 为卷积层数。图中对应的 f=2，K=4, 于是 R=16, 也就是当前时刻输出可以考虑到前面 16 个时刻步长输入。
与循环网络（RNN）相比，TCN 的缺点之一是它在推理期间占用更大的内存。计算下一个时间步需要整个原始序列。下图是一个经典的 TCN 模型结构（结合残差网络）：[链接4]
for i = 1:numBlocks
dilationFactor = 2^(i-1);
layers = [
convolution1dLayer(filterSize,numFilters,DilationFactor=dilationFactor,Padding="causal",Name="conv1_"+i)
layerNormalizationLayer
spatialDropoutLayer(dropoutFactor)
convolution1dLayer(filterSize,numFilters,DilationFactor=dilationFactor,Padding="causal")
layerNormalizationLayer
reluLayer
spatialDropoutLayer(dropoutFactor)
additionLayer(2,Name="add_"+i)];
% Add and connect layers.
lgraph =addLayers(lgraph,layers);
lgraph =connectLayers(lgraph,outputName,"conv1_"+i);
一般的 TCN 架构(如[1]中所述)由多个残差块组成，每个残差块包含两组具有相同扩张因子的扩张因果卷积层，然后是归一化层、ReLU 激活层和空间 dropout 层。
网络将每个块的输入加到该块的输出上(当输入和输出之间的通道数量不匹配时，对输入进行 1 × 1 的卷积)，并应用最终的激活函数。
循环神经网络: LSTM/Gru
循环网络的结构与前馈神经网络不同，它带有一个指向自身的环，用来表示它可以传递当前时刻处理的信息给下一时刻使用，我们选取 LSTM 来介绍，结构如下，
LSTM 之所以可以用于动态系统建模，是因为 LSTM 和动态系统有类似的特性：对于时刻 t, LSTM 使用当前网络的状态 (ht-1,ct-1) 和当前的输入 xt来计算网络输出 ht，同时更新当前网络的状态 (ht,ct)，ht 也叫输出状态，就是当前时刻 t 的 LSTM 网络的输出，ct 称为 cell state,包含学习来的历史时刻的状态信息。每个时间步 LSTM 都会对 ct进行更新：添加信息或移除信息，主要通过四个门函数(f,g,i,o)来实现，将前面的计算过程示意性的描述一下：
其中，上面方程中的 Wf,g,i,o, Rf,g,i,o, bf,g,i,o是这些门函数各自的可学习参数，主要包括针对输入 xt 的权重，针对上一时刻输出状态（同时也是当前时刻的输入）ht-1 的权重，以及偏置，这些可学习参数本身是无状态的，被所有时刻共享。训练的过程就是优化这些学习参数的过程。
总结一下：当前时刻的输出 ht不仅依赖于当前的输入xt，还依赖于当前的状态 (ht-1,ct-1)。它可以根据训练选择性的记住每一时刻的“重要”信息，并且利用这个信息作为状态，结合当前输入 xt 进行输出状态 ht预测。
电机温度预测示例
接下来我们使用数据结合 LSTM 模型来搭建永磁同步电机（PMSM）的一个代理模型（Surrogate Modeling），用于电机不同位置的温度预测。
类似 Demo [链接5]
理解数据集
数据集来自多片时长不同的数据。从数据中可以看到 PMSM 不同位置温度与电气系统和热系统对应的工况参数的相互影响。
同时环境温度的变化也会对电机不同位置温度以及相应需求扭矩所需的电流电压有影响。
我们使用环境温度、冷却液温度、电压、电流作为输入，输出为不同位置 PMSM 的温度。
数据预处理与特征工程
同样作为一个动态系统，我们通过对初始数据再进行不同尺度滑窗实现衍生特征生成，所有特征结合在一起作为 LSTM 的输入（尽管 LSTM 本身也具有考虑时间依赖关系的特性）。
% create derived features using raw voltages and currents derivedinputs =computedrivedfeatures(tt_data); % check the noise in the data tt_data=[tt_data derivedinputs]; Vnames=tt_data.Properties.VariableNames; s1=620;s2=2915;s3=4487;s4=8825; % preprocess exponentially weighted moving average [t_s1,t_s2,t_s3,t_s4]=preprocmovavg(tt_data,s1,s2,s3,s4,Vnames); % preprocess exponentially weighted moving variance [t_v1,t_v2,t_v3,t_v4]=preprocmovvar(tt_data,s1,s2,s3,s4,Vnames); % attach features to the original table predictors=[tt_data,t_s1,t_s2,t_s3,t_s4,t_v1,t_v2,t_v3,t_v4,tt_profileid]; responses=[tt(:,9:12) tt_profileid]; VResponsenames=responses.Properties.VariableNames;
准备训练数据和验证数据
holdOuts =[657258]; % define profiles that are withheld from training. [xtrain,ytrain]= prepareDataTrain(predictors,responses,holdOuts);
我们将使用第 58 条记录作为验证集，其中包括 4.64 小时的数据。
validationdata =58; [xvalidation, yvalidation]= prepareDataValidation(predictors,responses,validationdata); numResponses = size(ytrain{1},1); featureDimension = size(xtrain{1},1); numHiddenUnits=125;
上面的 DAG 网络可以对许多复杂的系统进行建模。通过上面的结构（左右两支）帮助我们对依赖于时间历史的物理行为动态行为以及含状态的行为进行建模。
查看预测结果
如图所示，红色和蓝色分别代表了实测数据和模型预测结果，右侧图像显示的是他们的残差，可以看到大部分误差在 1% 以内，预测效果比较理想。模型在瞬态变化较快和较慢的工况下都能和实测数据保持一致，说明模型也保持了一定的保真度。
将模型导出 Simulink
我们将训练好的模型保存为 .mat 文件，并将其导入 Simulink Deep Neural Network Predict 模块[链接9]，这样我们就有了一个只有 50Kbyte 大小可以预测温度的代理模型（Surrogate Model）用于仿真。
Neural Ordinary Differential Equations：神经网络 ODE
这是 NIPS 2018 年最佳论文[2]提出的一种新的网络层。当然这个要理论上介绍还是比较复杂，我们可以通过应用场景直观的解释一些 Neural ODE 如何实现动态系统建模。
试想我们有一个动态系统，因为系统动力学过于复杂，我们没有真正的物理模型，但我们可以不断地通过测量得到系统的初始状态 y(t0)，动态输 入 u(t0), u(t1),…, 与动态输出 y(t1),y(t2),…,y(tn)。接下来我们想是否可以实现这样一个微分方程：
使得这个微分方程正好代表了我们的系统，也就是说在系统输入 y(t0) 的初始条件下通过求解(例如使用 ode45)这个微分方程得出的解 yp(t) 和我们实测结果是吻合的。但如何基于已有的数据 y(t0),序列 u(t) 和序列 y(t) 得到 f 呢？如果我们将 f 用一个神经网络 F(θ) 替代，即
我们现在有数据 u(t), y(t)。我们不断地利用数据训练参数θ，使得上述方程的解 yp(t1), yp(t2),...,yp(tn) 与实测 y(t)是吻合的，那就可以得到这个动态系统的微分方程模型了，从而可以用于后续系统仿真与预测。如何理解 yp(t) 与实测 y(t) 的接近度,也就是损失函数? 我们简单介绍一下训练时 Loss 函数计算。对于方程，我们在知道系统初始状态 y(t0)，可以通过很多数值积分求解器（例如常用的 ode45）得到任何时刻的推断输出 yp(t) (当然前提是系统的输入 u(t)也是已知).
θ 是神经网络 F 的静态参数。对于所有时刻都是不变的。我们就可以方便得到损失函数的值
其中 L 可以是任何自定义的损失函数。于是我们可以进行参数 θ 的训练。当然关于梯度计算与反向传播会有相对复杂的数学推导，论文提出了伴随方法（Adjoint Method）来实现这一过程，此处不做详细论述。
对应于上述过程，MATLAB中提供了dlode45 [链接10]，用于建模方程右侧非线性函数F的同时，计算 ODE 的时序解。
即 dlode45 接收一个含参神经网络 F(θ)、需要计算输出结果的时刻序列 [t0, t1,…,tN]、系统的 t0 时刻的初始状态 yt0、神经网络参数的一组值，就可以计算出时刻 [t1,…,tN] 所对应的输出状态。
使用 Neural ODE 为系统建模示例
我们通过一个示例介绍如何使用Neural ODE为动态系统建模。[链接6]
我们就借用文章刚开始的简单二自由度线性系统，x' (t)=Ax(t), 其中A是一个 2x2 的矩阵。我们用这个已知的系统产生一些数据，利用这些数据来训练一个 Neural ODE 的方程，使得这个基于数据训练好的系统（Neural ODE 方程）能够接近已知的这个动态系统。
生成物理系统的数据 xTrain 作为真值
x0 = [2; 0];
A = [-0.1 -1; 1-0.1];
trueModel = @(t,y)A*y;
[~, xTrain] =ode45(trueModel, t, x0, odeOptions);
xTrain 两个自由度的可视化
定义和初始化神经网络 F(t,x(t),θ) 的参数 θ
neuralOdeParameters.fc1= struct;
sz = [hiddenSizestateSize];
neuralOdeParameters.fc1.Weights= initializeGlorot(sz, hiddenSize, stateSize);
neuralOdeParameters.fc1.Bias = initializeZeros([hiddenSize 1]);
neuralOdeParameters.fc1
ans = 包含以下字段的 struct:
Weights: [20×2 dlarray]
Bias: [20×1 dlarray]
同样
neuralOdeParameters.fc2
ans = 包含以下字段的 struct:
Weights: [2×20 dlarray]
Bias: [2×1 dlarray]
定义神经网络模型 F(t,x(t),θ) 函数
function y = odeModel(~,y,theta)
y =tanh(theta.fc1.Weights*y + theta.fc1.Bias);
y =theta.fc2.Weights*y + theta.fc2.Bias;
end
结合定义好的 F(t,x(t),θ) 作为 dlode45 的输入来构建代理模型函数
function X =model(tspan,X0,neuralOdeParameters)
X = dlode45(@odeModel,tspan,X0,neuralOdeParameters,DataFormat="CB");
end
定义模型梯度函数
主要用于训练过程计算损失以及对应待训练参数的梯度
function [gradients,loss] =modelGradients(tspan,dlX0,neuralOdeParameters,targets)
% Compute predictions.
dlX = model(tspan,dlX0,neuralOdeParameters);
% Compute L1 loss.
loss =l1loss(dlX,targets,NormalizationFactor="all-elements",DataFormat="CBT");
% Compute gradients.
gradients =dlgradient(loss,neuralOdeParameters);
end
训练模型
不断地迭代训练，创建 miniBatch,并进行损失函数计算和自动微分，通过调用 adam 求解器进行参数学习
for iter=1:numIter
% Create batch
[dlx0, targets] = createMiniBatch(numTrainingTimesteps,neuralOdeTimesteps, miniBatchSize, xTrain);
% Evaluatenetwork and compute gradients
[grads,loss] = dlfeval(@modelGradients,timesteps,dlx0,neuralOdeParameters,targets);
% Update network
[neuralOdeParameters,averageGrad,averageSqGrad] =adamupdate(neuralOdeParameters,grads,averageGrad,averageSqGrad,iter,...
learnRate,gradDecay,sqGradDecay);
% Plot loss
currentLoss =double(extractdata(loss));
测试模型
选取新的初始条件作为训练好的模型的输入，来进行和物理系统输出的对比
x0Pred1 =sqrt([2;2]);
x0Pred2 =[-1;-1.5];
x0Pred3 = [0;2];
x0Pred4 = [-2;0];
可以看到模型对于新的初始条件依然表现优异。因此神经网络 ODE 在构建动态系统上很有潜力，目前在发动机建模上也有一些示例应用。
NARX（nonlinear autoregressive network with exogenous inputs）反馈神经网络
在深度学习网络爆发之前，在浅层神经网络的场景中，NARX 反馈神经网络是经常用于动态系统建模的，它具有反馈连接，即输出 y(t) 依赖于系统之前时刻的输出 y(t-1)，y(t-2) 等等和输入。
其中网络 F 主要是通过一个前馈神经网络实现。
如下图，其中F即为第一层隐含层前馈网络，其中输入和输出可以是多维的。
跟其他反馈神经网络类似（前面介绍的 LSTM）,它的训练过程与推断过程有些区别。网络本身的输出需要被反馈到前馈神经网络的输入。在训练时，因为我们能够拿到整段输出真实的数据，因此，我们会用当前时刻真实的输出值作为训练时模型输入而不是反馈预测的输出，换句话说，在训练时我们会把网络作为开环去训练。
这有两个好处。一是前馈网络的输入更准确。第二，生成的网络具有纯前馈架构，静态反向传播变得可用。当我们进行未来多步推断时，因为这种情况我们只能用推断的数据进行下一时刻预测，所以这次我们才把网络闭环，用于推断。
磁悬浮系统的 NARX 建模示例
接下来我们使用 narx 神经网络来对动态系统进行建模的示例[链接7]。示例系统是一个磁悬浮系统。目标是控制悬浮在电磁铁上方的磁铁的位置，在电磁铁的位置上，磁铁受到限制，只能在垂直方向上移动，如下图所示，
系统的运动方程为：
其中 y(t) 是磁铁在电磁铁上方的距离，i(t) 是经过电磁铁的电流，M 是磁铁的质量，g 是重力常数。其中 β 为粘性摩擦系数，由磁体运动材料决定；α 为场强常数，由电磁铁上导线匝数和磁体强度决定。我们搜集了系统输入 u(t)-施加在电磁铁上的电压和系统输出 y(t)-永磁体的位置，对应两个时间序列。
搭建网络和准备数据
d1 = [1:2];
d2 = [1:2];
narx_net =narxnet(d1,d2,10); % 使用narxnet功能创建NARX开环串联网络，10个隐藏层神经元
[p,Pi,Ai,t] =preparets(narx_net,u,{},y); % 用preparets准备数据
其中 y(t) 是一个反馈信号，既是输入，也是输出，训练时我们既可以拿到当前时刻数据，也可以拿到后面时刻的数据，所以可以用于开环训练，当推断时我们会将输出接到输入作为闭环进行推断。
训练网络
narx_net =train(narx_net,p,t,Pi); % 训练网络得到训练好的网络
验证开环网络推断效果
yp =sim(narx_net,p,Pi);
e =cell2mat(yp)-cell2mat(t);
可以看到误差很小。因为我们用的开环训练，所以推断结果是用前面时刻的真实输出数据（而非推断输出反馈），所以这里的误差是 one-step-ahead 推断误差。
测试闭环推断效果
如果要看网络真实准确度的表现，需要将开环的输出作为反馈接到输入，然后进行多步预测。
narx_net_closed =closeloop(narx_net);
现在可以使用闭环执行 900 个时间步的迭代预测。在这个闭环网络中只需要两个初始输入和两个初始输出作为初始条件。
y1 = y(1700:2600);
u1 = u(1700:2600);
[p1,Pi1,Ai1,t1] = preparets(narx_net_closed,u1,{},y1);
yp1 =narx_net_closed(p1,Pi1,Ai1);
从闭环预测的结果看，蓝线是磁铁的实际位置，红线是 NARX 神经网络预测的位置。即使网络预测的时间步预测了 900 步，预测依然是非常准确的。闭环多步预测准确的前提就是开环下单步推断误差要小。
总结
本文主要介绍了动态系统的特性和用于动态系统建模的神经网络模型从而可以实现模型降阶（Reduced Order Modeling），包括前馈神经网络，TCN，循环神经网络，神经网络 ODE, NARX 网络和相应的一些手段，并结合示例与场景进行了说明这些手段的有效性。
在后续文章我会来介绍系统辨识（System Identification）的一些示例,欢迎继续关注。
本文中涉及到的多个 demo，数据以及脚本文件，若您感兴趣，可以通过扫描填写下面这个反馈问卷，或点击”阅读原文“进一步获取这些链接。
获取文中示例链接
参考文献
[1] Vidal, C., Kollmeyer, P., Naguib, M., Malysz, P. et al., “Robust xEV Battery State-of-Charge Estimator Design Using a Feedforward Deep Neural Network,” SAE Technical Paper 2020-01-1181, 2020, doi:10.4271/2020-01-1181.
[2] Ricky T. Q. Chen*, Yulia Rubanova*, Jesse Bettencourt*, David Duvenaud University of Toronto, Vector Institute “Neural Ordinary Differential Equations”
编辑：谢雅洁 校对 ：向映姣