我們將訓練好的模型保存為 .mat 文件，并將其導入 Simulink Deep Neural Network Predict 模塊[鏈接9]，這樣我們就有了一個只有 50Kbyte 大小可以預測溫度的代理模型（Surrogate Model）用于仿真。

Neural Ordinary Differential Equations：神經(jīng)網(wǎng)絡 ODE

這是 NIPS 2018 年最佳論文[2]提出的一種新的網(wǎng)絡層。當然這個要理論上介紹還是比較復雜，我們可以通過應用場景直觀的解釋一些 Neural ODE 如何實現(xiàn)動態(tài)系統(tǒng)建模。

試想我們有一個動態(tài)系統(tǒng)，因為系統(tǒng)動力學過于復雜，我們沒有真正的物理模型，但我們可以不斷地通過測量得到系統(tǒng)的初始狀態(tài) y(t0)，動態(tài)輸 入 u(t0), u(t1),…, 與動態(tài)輸出  y(t1),y(t2),…,y(tn)。接下來我們想是否可以實現(xiàn)這樣一個微分方程：

使得這個微分方程正好代表了我們的系統(tǒng)，也就是說在系統(tǒng)輸入 y(t0) 的初始條件下通過求解(例如使用 ode45)這個微分方程得出的解 yp(t) 和我們實測結(jié)果是吻合的。但如何基于已有的數(shù)據(jù) y(t0),序列 u(t) 和序列 y(t) 得到 f 呢？如果我們將 f 用一個神經(jīng)網(wǎng)絡 F(θ) 替代，即

我們現(xiàn)在有數(shù)據(jù) u(t), y(t)。我們不斷地利用數(shù)據(jù)訓練參數(shù)θ，使得上述方程的解 yp(t1), yp(t2),...,yp(tn) 與實測 y(t)是吻合的，那就可以得到這個動態(tài)系統(tǒng)的微分方程模型了，從而可以用于后續(xù)系統(tǒng)仿真與預測。如何理解  yp(t)  與實測 y(t) 的接近度,也就是損失函數(shù)? 我們簡單介紹一下訓練時 Loss 函數(shù)計算。對于方程，我們在知道系統(tǒng)初始狀態(tài) y(t0)，可以通過很多數(shù)值積分求解器（例如常用的 ode45）得到任何時刻的推斷輸出  yp(t) (當然前提是系統(tǒng)的輸入 u(t)也是已知).

即

θ 是神經(jīng)網(wǎng)絡 F 的靜態(tài)參數(shù)。對于所有時刻都是不變的。我們就可以方便得到損失函數(shù)的值

其中 L 可以是任何自定義的損失函數(shù)。于是我們可以進行參數(shù) θ 的訓練。當然關于梯度計算與反向傳播會有相對復雜的數(shù)學推導，論文提出了伴隨方法（Adjoint Method）來實現(xiàn)這一過程，此處不做詳細論述。

對應于上述過程，MATLAB中提供了dlode45 [鏈接10]，用于建模方程右側(cè)非線性函數(shù)F的同時，計算 ODE 的時序解。

即 dlode45 接收一個含參神經(jīng)網(wǎng)絡 F(θ)、需要計算輸出結(jié)果的時刻序列 [t0, t1,…,tN]、系統(tǒng)的 t0 時刻的初始狀態(tài) yt0、神經(jīng)網(wǎng)絡參數(shù)的一組值，就可以計算出時刻 [t1,…,tN] 所對應的輸出狀態(tài)。

使用 Neural ODE 為系統(tǒng)建模示例

我們通過一個示例介紹如何使用Neural ODE為動態(tài)系統(tǒng)建模。[鏈接6]

我們就借用文章剛開始的簡單二自由度線性系統(tǒng)，x' (t)=Ax(t), 其中A是一個 2x2 的矩陣。我們用這個已知的系統(tǒng)產(chǎn)生一些數(shù)據(jù)，利用這些數(shù)據(jù)來訓練一個 Neural ODE 的方程，使得這個基于數(shù)據(jù)訓練好的系統(tǒng)（Neural ODE 方程）能夠接近已知的這個動態(tài)系統(tǒng)。

生成物理系統(tǒng)的數(shù)據(jù) xTrain 作為真值

x0 = [2; 0];

A = [-0.1 -1; 1-0.1];

trueModel = @(t,y)A*y;

[~, xTrain] =ode45(trueModel, t, x0, odeOptions);

xTrain 兩個自由度的可視化

定義和初始化神經(jīng)網(wǎng)絡 F(t,x(t),θ) 的參數(shù) θ

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]

定義神經(jīng)網(wǎng)絡模型 F(t,x(t),θ) 函數(shù)

function y = odeModel(~,y,theta)

    y =tanh(theta.fc1.Weights*y + theta.fc1.Bias);

    y =theta.fc2.Weights*y + theta.fc2.Bias;

end

結(jié)合定義好的 F(t,x(t),θ) 作為 dlode45 的輸入來構(gòu)建代理模型函數(shù)

function X =model(tspan,X0,neuralOdeParameters)

    X = dlode45(@odeModel,tspan,X0,neuralOdeParameters,DataFormat="CB");

end

定義模型梯度函數(shù)

主要用于訓練過程計算損失以及對應待訓練參數(shù)的梯度

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

訓練模型

不斷地迭代訓練，創(chuàng)建 miniBatch,并進行損失函數(shù)計算和自動微分，通過調(diào)用 adam 求解器進行參數(shù)學習

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));

…

測試模型

選取新的初始條件作為訓練好的模型的輸入，來進行和物理系統(tǒng)輸出的對比

x0Pred1 =sqrt([2;2]);

x0Pred2 =[-1;-1.5];

x0Pred3 = [0;2];

x0Pred4 = [-2;0];

可以看到模型對于新的初始條件依然表現(xiàn)優(yōu)異。因此神經(jīng)網(wǎng)絡 ODE 在構(gòu)建動態(tài)系統(tǒng)上很有潛力，目前在發(fā)動機建模上也有一些示例應用。