使用LabVIEW調(diào)用IR模型解薛定諤方程

作者：王立奇英特爾邊緣計算創(chuàng)新大使

一、PINN—加入物理約束的神經(jīng)網(wǎng)絡(luò)

基于物理信息的神經(jīng)網(wǎng)絡(luò)（Physics-informed Neural Network，簡稱PINN），是一類用于解決有監(jiān)督學(xué)習(xí)任務(wù)的神經(jīng)網(wǎng)絡(luò)，它不僅能夠像傳統(tǒng)神經(jīng)網(wǎng)絡(luò)一樣學(xué)習(xí)到訓(xùn)練數(shù)據(jù)樣本的分布規(guī)律，而且能夠?qū)W習(xí)到數(shù)學(xué)方程描述的物理定律。與純數(shù)據(jù)驅(qū)動的神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)相比，PINN在訓(xùn)練過程中施加了物理信息約束，因而能用更少的數(shù)據(jù)樣本學(xué)習(xí)到更具泛化能力的模型。本文主要解析這種神經(jīng)網(wǎng)絡(luò)以及相關(guān)應(yīng)用。

2019年，來自布朗大學(xué)應(yīng)用數(shù)學(xué)的研究團隊提出了一種用物理方程作為運算限制的“物理激發(fā)的神經(jīng)網(wǎng)絡(luò)”(PINN)并發(fā)表在了計算物理學(xué)領(lǐng)域權(quán)威雜志《計算物理學(xué)期刊》(Journal of Computational Physics)上。這篇論文一經(jīng)發(fā)表就獲得了大量關(guān)注。這篇論文因為代碼體系的完整性使得開發(fā)人員們很容易上手把相關(guān)的學(xué)習(xí)框架應(yīng)用到不同領(lǐng)域上去。所以在發(fā)表不久之后，一系列不同的PINN也被其他研究者開發(fā)出來。甚至可以不夸張的說，PINN是目前AI物理領(lǐng)域論文中最常見到的框架和詞匯之一。

算法描述

而所謂的物理神經(jīng)網(wǎng)絡(luò)，其實就是把物理方程作為限制加入神經(jīng)網(wǎng)絡(luò)中使訓(xùn)練的結(jié)果滿足物理規(guī)律。而這個所謂的限制是怎么實現(xiàn)的？其實就是通過把物理方程的迭代前后的差值加到神經(jīng)網(wǎng)絡(luò)的損失函數(shù)里面去，讓物理方程也“參與”到了訓(xùn)練過程。這樣，神經(jīng)網(wǎng)絡(luò)在訓(xùn)練迭代時候優(yōu)化的不僅僅是網(wǎng)絡(luò)自己的損失函數(shù)，還有物理方程每次迭代的差，使得最后訓(xùn)練出來的結(jié)果就滿足物理規(guī)律了。

二、DeepXDE

DeepXDE 由 Lu Lu 在布朗大學(xué) George Karniadakis 教授的指導(dǎo)下于 2018 年夏季至 2020年夏季開發(fā)，并得到 PhILM 的支持。DeepXDE 最初是在布朗大學(xué)的 Subversion 中自行托管的，名稱為 SciCoNet（科學(xué)計算神經(jīng)網(wǎng)絡(luò)）。2019 年 2 月 7 日，SciCoNet 從 Subversion 遷移到 GitHub，更名為 DeepXDE。

DeepXDE特性

DeepXDE 已經(jīng)實現(xiàn)了如上所示的許多算法，并支持許多特性：

復(fù)雜的域幾何圖形，沒有專制網(wǎng)格生成。原始幾何形狀是間隔、三角形、矩形、多邊形、圓盤、長方體和球體。其他幾何可以使用三個布爾運算構(gòu)建為構(gòu)造實體幾何(CSG)：并集、差集和交集。

多物理場，即（時間相關(guān)的）耦合偏微分方程。

5 種類型的邊界條件(BC)：Dirichlet、Neumann、Robin、周期性和一般 BC，可以在任意域或點集上定義。

不同的神經(jīng)網(wǎng)絡(luò)，例如（堆疊/非堆疊）全連接神經(jīng)網(wǎng)絡(luò)、殘差神經(jīng)網(wǎng)絡(luò)和（時空）多尺度傅里葉特征網(wǎng)絡(luò)。

6種抽樣方法：均勻抽樣、偽隨機抽樣、拉丁超立方抽樣、Halton序列、Hammersley序列、Sobol序列。

訓(xùn)練點可以在訓(xùn)練期間保持不變，也可以每隔一定的迭代重新采樣一次。

方便保存訓(xùn)練期間的模型，并加載訓(xùn)練好的模型。

使用 dropout 的不確定性量化。

許多不同的（加權(quán)）損失、優(yōu)化器、學(xué)習(xí)率計劃、指標(biāo)等回調(diào)，用于在訓(xùn)練期間監(jiān)控模型的內(nèi)部狀態(tài)和統(tǒng)計信息，例如提前停止。

使用戶代碼緊湊，與數(shù)學(xué)公式非常相似。

DeepXDE 的所有組件都是松耦合的，因此 DeepXDE 結(jié)構(gòu)良好且高度可配置。

可以輕松自定義 DeepXDE 以滿足新的需求。

數(shù)值算例

（1）問題設(shè)置

我們將求解由下式給出的非線性薛定諤方程：

周期性邊界條件為：

初始條件為：

Deepxde 只使用實數(shù)，因此我們需要明確拆分復(fù)數(shù) PDE 的實部和虛部。

代替單個殘差：

我們得到兩個（實值）殘差：

其中 u ( x , t )和 v ( x , t )分別表示h的實部和虛部。

import numpy as np
import deepxde as dde

#用于繪圖

import matplotlib.pyplot as plt
from scipy.interpolate import griddata


x_lower = -5
x_upper = 5
t_lower = 0
t_upper = np.pi / 2

左滑查看更多

#創(chuàng)建2D域（用于繪圖和輸入）

x = np.linspace(x_lower, x_upper, 256)
t = np.linspace(t_lower, t_upper, 201)
X, T = np.meshgrid(x, t)

左滑查看更多

#整個域變平

X_star = np.hstack((X.flatten()[:, None], T.flatten()[:, None]))

左滑查看更多

#空間和時間域/幾何（對于deepxde模型）

space_domain = dde.geometry.Interval(x_lower, x_upper)
time_domain = dde.geometry.TimeDomain(t_lower, t_upper)
geomtime = dde.geometry.GeometryXTime(space_domain, time_domain)

左滑查看更多

#損失的“物理信息”部分

def pde(x, y):
  """
  INPUTS:
    x: x[:,0] 是 x 坐標(biāo)
      x[:,1] 是 t 坐標(biāo)
    y: 網(wǎng)絡(luò)輸出，在這種情況下:
      y[:,0] 是 u(x,t) 實部
      y[:,1] 是 v(x,t) 虛部
  OUTPUT:
    標(biāo)準(zhǔn)形式的 pde，即必須為零的東西
  """


  u = y[:, 0:1]
  v = y[:, 1:2]

左滑查看更多

#在'jacobian'中，i是輸出分量，j是輸入分量

 u_t = dde.grad.jacobian(y, x, i=0, j=1)
  v_t = dde.grad.jacobian(y, x, i=1, j=1)


  u_x = dde.grad.jacobian(y, x, i=0, j=0)
  v_x = dde.grad.jacobian(y, x, i=1, j=0)

左滑查看更多

#在“hessian”中，i 和 j 都是輸入分量。（Hessian 原則上可以是 d^2y/dxdt、d^2y/d^2x 等）

#輸出組件由“組件”選擇

 u_xx = dde.grad.hessian(y, x, component=0, i=0, j=0)
  v_xx = dde.grad.hessian(y, x, component=1, i=0, j=0)


  f_u = u_t + 0.5 * v_xx + (u ** 2 + v ** 2) * v
  f_v = v_t - 0.5 * u_xx - (u ** 2 + v ** 2) * u


  return [f_u, f_v]

左滑查看更多

#邊界條件和初始條件

#周期性邊界條件

bc_u_0 = dde.PeriodicBC(
  geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=0, component=0
)
bc_u_1 = dde.PeriodicBC(
  geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=1, component=0
)
bc_v_0 = dde.PeriodicBC(
  geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=0, component=1
)
bc_v_1 = dde.PeriodicBC(
  geomtime, 0, lambda _, on_boundary: on_boundary, derivative_order=1, component=1
)

左滑查看更多

#初始條件

def init_cond_u(x):
  "2 sech(x)"
  return 2 / np.cosh(x[:, 0:1])
  
def init_cond_v(x):
  return 0




ic_u = dde.IC(geomtime, init_cond_u, lambda _, on_initial: on_initial, component=0)
ic_v = dde.IC(geomtime, init_cond_v, lambda _, on_initial: on_initial, component=1)
data = dde.data.TimePDE(
  geomtime,
  pde,
  [bc_u_0, bc_u_1, bc_v_0, bc_v_1, ic_u, ic_v],
  num_domain=10000,
  num_boundary=20,
  num_initial=200,
  train_distribution="pseudo",
)

左滑查看更多

#網(wǎng)絡(luò)架構(gòu)

net = dde.maps.FNN([2] + [100] * 4 + [2], "tanh", "Glorot normal")


model = dde.Model(data, net)

左滑查看更多

Adam優(yōu)化

#強烈建議使用GPU加速系統(tǒng).

model.compile("adam", lr=1e-3, loss="MSE")
model.train(epochs=1000, display_every=100)

左滑查看更多

L-BFGS優(yōu)化

dde.optimizers.config.set_LBFGS_options(
  maxcor=50,
  ftol=1.0 * np.finfo(float).eps,
  gtol=1e-08,
  maxiter=1000,
  maxfun=1000,
  maxls=50,
)
model.compile("L-BFGS")
model.train()
Compiling model...
'compile' took 0.554160 s


Training model...


Step   Train loss                                     Test loss                                      Test metric
1000   [1.25e-02, 1.32e-02, 2.72e-06, 3.28e-04, 5.27e-06, 4.65e-06, 2.24e-02, 1.31e-03]  [1.25e-02, 1.32e-02, 2.72e-06, 3.28e-04, 5.27e-06, 4.65e-06, 2.24e-02, 1.31e-03]  [] 
2000   [7.03e-04, 7.62e-04, 6.76e-06, 1.33e-05, 2.88e-07, 8.49e-06, 4.01e-04, 3.86e-05]                                              
INFOOptimization terminated with:
 Message: STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT
 Objective function value: 0.001928
 Number of iterations: 945
 Number of functions evaluations: 1001
2001   [7.18e-04, 7.43e-04, 6.27e-06, 1.23e-05, 2.94e-07, 8.89e-06, 4.01e-04, 3.82e-05]  [7.18e-04, 7.43e-04, 6.27e-06, 1.23e-05, 2.94e-07, 8.89e-06, 4.01e-04, 3.82e-05]  [] 


Best model at step 2001:
 train loss: 1.93e-03
 test loss: 1.93e-03
 test metric: []


'train' took 179.449384 s












(,
 )

左滑查看更多

#做預(yù)測

prediction = model.predict(X_star, operator=None)


u = griddata(X_star, prediction[:, 0], (X, T), method="cubic")
v = griddata(X_star, prediction[:, 1], (X, T), method="cubic")


h = np.sqrt(u ** 2 + v ** 2)

左滑查看更多

#繪制預(yù)測

fig, ax = plt.subplots(3)


ax[0].set_title("Results")
ax[0].set_ylabel("Real part")
ax[0].imshow(
  u.T,
  interpolation="nearest",
  cmap="viridis",
  extent=[t_lower, t_upper, x_lower, x_upper],
  origin="lower",
  aspect="auto",
)
ax[1].set_ylabel("Imaginary part")
ax[1].imshow(
  v.T,
  interpolation="nearest",
  cmap="viridis",
  extent=[t_lower, t_upper, x_lower, x_upper],
  origin="lower",
  aspect="auto",
)
ax[2].set_ylabel("Amplitude")
ax[2].imshow(
  h.T,
  interpolation="nearest",
  cmap="viridis",
  extent=[t_lower, t_upper, x_lower, x_upper],
  origin="lower",
  aspect="auto",
)


plt.show()

左滑查看更多

三、生成OpenVINO IR模型

如需使用LabVIEW OpenVINO推理deepxde的模型，需要做以下幾步：

將deepxde模型轉(zhuǎn)成onnx模型作為中間件：

代碼中，我們首先要創(chuàng)建一個和X_star的形狀相同的X2作為參考輸入（X2必須是cpu上的矩陣），其次是新建一個基于cpu的mode2，其權(quán)重和model完全一致。最后使用torch.onnx.export生成onnx模型。

使用命令行將onnx模型轉(zhuǎn)為IR模型：

mo --input_model Schrodinger.onnx --input_shape "[256,2]"

左滑查看更多

命令中，我們將輸入的形狀設(shè)置為256*2，其中第0列為位置x，第1列為時間t。256為我們設(shè)置的每個時間點需要推理的x的點數(shù)，可以為任意長度。

完成上述兩個步驟后，我們可以看到文件夾里多了三個文件Schrodinger.onnx、Schrodinger.xml和Schrodinger.bin。其中onnx為pytorch生成的onnx模型，xml和bin為mo生成的OpenVINO IR模型。

四、LabVIEW調(diào)用IR模型

LabVIEW是NI推出的圖形化編程環(huán)境，在科研、工業(yè)測控領(lǐng)域有著廣泛的應(yīng)用。使用LabVIEW的波形圖、強度圖等控件，可以使用戶更直觀的觀測自己模型的訓(xùn)練結(jié)果。

使用我們編寫好的Schrodinger_OpenVINO.vi，即可快速推理Schrodinger方程的OpenVINO模型。以下是使用LabVIEW調(diào)用模型的步驟：

1. 初始化模型：使用LoadIR.vi調(diào)用xml和bin文件，推理引擎可使用CPU或GPU（Intel核顯或獨顯）；

2. 初始化位置。根據(jù)我們訓(xùn)練時用的參數(shù)，需要模擬-5~5范圍內(nèi)任意數(shù)量的位置點（本案例中使用256個位置點）。

3. 初始化時間，根據(jù)我們訓(xùn)練的參數(shù)，需要模擬0~pi/2范圍內(nèi)任意數(shù)量的時間點（本案例中使用201個時間點。

4. 推理模型時。每次循環(huán)將某一個時間點復(fù)制256份，和所有的256個位置點組成2*256的二維數(shù)組，然后轉(zhuǎn)置成256*2的數(shù)組，輸入至網(wǎng)絡(luò)中。

5. 獲取輸出。輸出為256*2，但和輸入的意義不同：其中第0列為波函數(shù)h的實部，第1列為h的虛部。因此需要將結(jié)果轉(zhuǎn)置成2*256，后，方能使用波形圖（Waveform Graph）顯示兩條曲線。

6. 獲取整個時間段、所有位置的結(jié)果。使用For循環(huán)索引，保存所有時間段實部、虛部和幅值的結(jié)果，并顯示在強度圖上。

7. 使用release.vi關(guān)閉OpenVINO模型。

運行程序。我們使用Intel A770獨立顯卡(設(shè)備名稱為GPU.1)運行該模型。運行結(jié)果如下圖：

左邊的波形圖為每一個時間點的波函數(shù)實部和虛部的曲線，右邊的強度圖為所有時間點的實部、虛部和幅值分布。結(jié)果表明，使用OpenVINO推理物理信息神經(jīng)網(wǎng)絡(luò)模型，和deepxde的結(jié)果完全一致。此外，在A770上模型也得到了最大加速，并節(jié)省了大量的cpu資源。

公司簡介

上海儀酷智能科技有限公司是國內(nèi)領(lǐng)先的人工智能軟件提供商和機器人應(yīng)用開發(fā)商，開發(fā)了低代碼可重構(gòu)的人工智能定制平臺，包括基于邊緣計算的高性能、高性價比工業(yè)硬件，以及跨平臺的圖形化深度學(xué)習(xí)視覺軟件，滿足智能制造、智慧醫(yī)療、智慧教育等多個領(lǐng)域的定制化需求；同時為廣大院校用戶搭建了人工智能教育軟硬件平臺（語音交互套件、無人駕駛套件、機械臂套件、以及支持所有硬件的云端實驗平臺）；公司是百度、軟銀機器人、NI、Intel的戰(zhàn)略合作伙伴，科大訊飛創(chuàng)投生態(tài)圈下企業(yè)，連續(xù)多年上海交通大學(xué)優(yōu)秀合作伙伴，上海市高新技術(shù)企業(yè)，姑蘇領(lǐng)軍人才,且入駐百度飛槳人工智能產(chǎn)業(yè)賦能中心，擁有70余項專利及著作權(quán)，總部位于上海，于蘇州、安慶設(shè)立分公司。

本文使用LabVIEW調(diào)用IR模型解薛定諤方程，主要使用了儀酷智能LabVIEW圖形化AI視覺工具包和LabVIEW圖形化OpenVINO工具包，大家可以在儀酷智能公司官網(wǎng)下載相應(yīng)軟件工具包。