Physics-informed Neural Networks (PINN)
http://iai.postech.ac.kr/
Industrial AI Lab at POSTECH
Table of Contents
1. Why Deep Learning Needs Physics?Ā¶
Why do data-driven āblack-boxā methods fail?
- May output result that is physically inconsistent
- Easy to learn spurious relationships that look good only on training and test data
- Can lead to poor generalization outside the available data (out-of-sample prediction tasks)
- Interpretability is absent
- Discovering the mechanism of an underlying process is crucial for scientific advancements
Physics-Informed Neural Networks (PINNs)
- Take full advantage of data science methods with the accumulated prior knowledge of scientific theories $\rightarrow$ Improve predictive performance
- Integration of domain knowledge to overcome the issue of imbalanced data & data shortage
1.1. Taxonomy of Informed Deep LearningĀ¶
1.2. Multilayer Feedforward Networks are Universal ApproximatorsĀ¶
- The Universal Approximation Theorem
- Neural Networks are capable of approximating any Borel measurable function
- Neural Networks (1989)
1.3. Neural Networks for Solving Differential EquationsĀ¶
- Neural Algorithm for Solving Differential Equations
- Journal of Computational Physics (1990)
- Neural minimization for finite difference equation
- ANN for ODE and PDE
- IEEE on Neural Networks (1998)
1.4. Journal of Computational Physics (2019)Ā¶
1.5. Nature Reviews Physics (2021)Ā¶
2. Architecture of Physics-informed Neural Networks (PINN)Ā¶
- NN as an universal function approximator
- Given
- Some measured data from initial and boundary conditions
- ODE or PDE
$$
\frac{\partial u}{\partial t} - \lambda \frac{\partial^2 u}{\partial x^2} = 0
$$
- Adding constraints for regularization
- Regularized by physics, but matched with sparse data
$$
\begin{align*}
\mathcal{L} &= \omega_{\text{data}} \mathcal{L}_{\text{data}} + \omega_{\text{PDE}} \mathcal{L}_{\text{PDE}}\\\\\\
\mathcal{L}_{\text{data}} &= \frac{1}{N_{\text{data}}} \sum \left(u(x,t) - u \right)^2 \\
\mathcal{L}_{\text{PDE}} &= \frac{1}{N_{\text{PDE}}} \sum \left(\frac{\partial u}{\partial t} - \lambda \frac{\partial^2 u}{\partial x^2} \right)^2
\end{align*}
$$
3. Methond for Solving ODE with Neural NetworksĀ¶
3.1. BackgroundĀ¶
This is a result first due to Lagaris et.al from 1998. The idea is to solve differential equations using neural networks by representing the solution by a neural network and training the resulting network to satisfy the conditions required by the differential equation.
Consider a system of ordinary differential equations
$$ {u^\prime} = f(u,t) $$with $t \in [0, 1]$ and a known initial condition $u(0) = u_0$.
To solve this, we approximate the solution by a neural network:
$$ {\text{NN}(t)} \approx {u(t)} $$If $\text{NN}(t)$ was the true solution, then it would hold that $\text{NN}'(t) = f\left(\text{NN}(t),t \right)$ for all $t$. Thus we turn this condtion into our loss function. This motivates the loss function.
$${L(\omega)} = \sum_{i} \left(\frac{d \text{NN}(t_i)}{dt}-f\left(\text{NN}(t_i),t_i \right)\right)^2$$
The choice of $t_i$ could be done in many ways: it can be random, it can be a grid, etc. Anyways, when this loss function is minimized (gradients computed with standard reverse-mode automatic differentiation), then we have that $\frac{dNN(t_i)}{dt} \approx f(\text{NN}(t_i),t_i)$ and thus $\text{NN}(t)$ approximately solves the differential equation.
Note that we still have to handle the initial condition. One simple way to do this is to add an initial condition term to the cost function.
$${L(\omega)} = \sum_{i} \left(\frac{d \text{NN}(t_i)}{dt}-f(\text{NN}(t_i),t_i)\right)^2 + (\text{NN}(0)-u_0)^2 $$
where $\omega$ are the parameters that define the neural network $\text{NN}$ that approximates $u$. Thus this reduces down, once again, to simply finding weights which minimize a loss function !
3.2. Lab 1: Simple ExampleĀ¶
- Let's look at the ODE:
$$\frac{du}{dt} = \cos2\pi t$$
- Initial condition:
$$u(0) = 1$$
- The exact solution:
$$u(t) = \frac{1}{2\pi}\sin2\pi t + 1$$
- Make a neural network and loss functions like below:
3.2.0. SetupĀ¶
- Using PyTorch
- Install and Setup
- TensorFlow >= 2.2.0 and TensorFlow Probability >= 0.10.0
- If you have TensorFlow 2.2.0, TensorFlow Probability == 0.10.0 should be installed
- (TensorFlow 2.3.0 matches with TensorFlow Probability 0.11.0 and so on)
3.2.1. Import LibraryĀ¶
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
import random
3.2.2. Define Network and Hyper-parametersĀ¶
NN = tf.keras.models.Sequential([
tf.keras.layers.Input((1,)),
tf.keras.layers.Dense(units = 32, activation = 'tanh'),
tf.keras.layers.Dense(units = 32, activation = 'tanh'),
tf.keras.layers.Dense(units = 32, activation = 'tanh'),
tf.keras.layers.Dense(units = 1)
])
NN.summary()
optm = tf.keras.optimizers.Adam(learning_rate = 0.001)
3.2.3. Define ODE SystemĀ¶
- ODE loss:
$$L_{\text{ODE}} = \frac{1}{n}\sum_{i=1}^{n} \left(\frac{d \text{NN}(t_i)}{dt}- \cos2\pi t_i\right)^2$$
- Initial condition loss:
$$L_{IC} = \frac{1}{n}\sum_{i=1}^{n} \left({\text{NN}(0)}- 1\right)^2$$
- Total loss:
$$L_{\text{Total}} = L_{\text{ODE}} + L_{\text{IC}}$$
def ode_system(t, net):
t = t.reshape(-1,1)
t = tf.constant(t, dtype = tf.float32)
t_0 = tf.zeros((1,1))
one = tf.ones((1,1))
with tf.GradientTape() as tape:
tape.watch(t)
u = net(t)
u_t = tape.gradient(u, t)
ode_loss = u_t - tf.math.cos(2*np.pi*t)
IC_loss = net(t_0) - one
square_loss = tf.square(ode_loss) + tf.square(IC_loss)
total_loss = tf.reduce_mean(square_loss)
return total_loss
3.2.4. Train & PredictionĀ¶
train_t = (np.array([0., 0.025, 0.475, 0.5, 0.525, 0.9, 0.95, 1., 1.05, 1.1, 1.4, 1.45, 1.5, 1.55, 1.6, 1.95, 2.])).reshape(-1, 1)
train_loss_record = []
for itr in range(6000):
with tf.GradientTape() as tape:
train_loss = ode_system(train_t, NN)
train_loss_record.append(train_loss)
grad_w = tape.gradient(train_loss, NN.trainable_variables)
optm.apply_gradients(zip(grad_w, NN.trainable_variables))
if itr % 1000 == 0:
print(train_loss.numpy())
plt.figure(figsize = (10,8))
plt.plot(train_loss_record)
plt.show()
test_t = np.linspace(0, 2, 100)
train_u = np.sin(2*np.pi*train_t)/(2*np.pi) + 1
true_u = np.sin(2*np.pi*test_t)/(2*np.pi) + 1
pred_u = NN.predict(test_t).ravel()
plt.figure(figsize = (10,8))
plt.plot(train_t, train_u, 'ok', label = 'Train')
plt.plot(test_t, true_u, '-k',label = 'True')
plt.plot(test_t, pred_u, '--r', label = 'Prediction')
plt.legend(fontsize = 15)
plt.xlabel('t', fontsize = 15)
plt.ylabel('u', fontsize = 15)
plt.show()
!pip install deepxde
- Change DeepXDE backends
- DeepXDE supports backends for TensorFlow 1.x, TensorFlow 2.x, and PyTorch. We will change it into Pytorch
# from deepxde.backend.set_default_backend import set_default_backend
# set_default_backend("tensorflow")
4.2. Problem Setup AgainĀ¶
- ODE
$$\frac{du}{dt} = \cos2\pi t$$
- Boundary condition:
$$u(0) = 1$$
- The exact solution:
$$u(t) = \frac{1}{2\pi}\sin2\pi t + 1$$
4.2.1. Import LibraryĀ¶
import tensorflow as tf
import deepxde as dde
import numpy as np
import matplotlib.pyplot as plt
import math as m
4.2.2. Define ODE SystemĀ¶
pi = tf.constant(m.pi)
def ode_system(t, u):
du_t = dde.grad.jacobian(u, t)
return du_t - tf.math.cos(2*pi*t)
4.2.3. Define Initial ConditionĀ¶
def boundary(t, on_initial):
return on_initial and np.isclose(t[0], 0)
4.2.4. Define Geometry & Implement Initial ConditionĀ¶
geom = dde.geometry.TimeDomain(0, 2)
ic = dde.IC(geom, lambda t: 1, boundary)
# Reference solution to compute the error
def true_solution(t):
return np.sin(2*np.pi*t)/(2*np.pi) + 1
data = dde.data.PDE(geom,
ode_system,
ic,
num_domain = 30,
num_boundary = 2,
solution = true_solution,
num_test = 100)
4.2.5. Define Network and Hyper-parametersĀ¶
layer_size = [1] + [32] + [1]
activation = "tanh"
initializer = "Glorot uniform"
NN = dde.maps.FNN(layer_size, activation, initializer)
model = dde.Model(data, NN)
model.compile("adam", lr = 0.001)
4.2.6. Train & PredictionĀ¶
losshistory, train_state = model.train(epochs = 6000)
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
5. Lab 3: Spring and Mass SystemĀ¶
5.1. Problem SetupĀ¶
- We will solve a spring coupled masses problem:
- Problem properties
- Initial properties
- Partial differential equations
$$m_1 \ddot{x_1} + k_1 (x_1 - L_1) - k_2 (x_2 -x_1 - L_2) = 0$$
$$m_2 \ddot{x_2} + k_2 (x_2 -x_1 - L_2) = 0$$
- One Dirichlet initial condition:
$$x_1(0) = 0.5, \quad x_2(0) = 3.25$$
- One Neumann initial condition:
$$\dot{x_1}(0) = 0, \quad \dot{x_2}(0) = 0$$
- Make a neural network and loss functions like below :
5.2. Solve the Spring Coupled Masses ProblemĀ¶
5.2.1. Import LibraryĀ¶
import deepxde as dde
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
5.2.2. Initial ConstantsmĀ¶
# Masses
m1 = 1.0
m2 = 1.0
# Spring constants
k1 = 5.0
k2 = 2.0
# Natural lengths
L1 = 0.5
L2 = 0.5
# Initial conditions
# d1_0 and d2_0 are the initial displacements; v1_0 and v2_0 are the initial velocities
d1_0 = 0.5
d2_0 = 3.25
v1_0 = 0.0
v2_0 = 0.0
b1 = 0
b2 = 0
# maximum time to simulate
t_max = 5
5.2.3. Define PDE SystemĀ¶
def dy(t, x):
return dde.grad.jacobian(x, t)
def pde(t, x):
# mass 1 location
x_1 = x[:, [0]]
# mass 2 location
x_2 = x[:, [1]]
dx1_tt = dde.grad.hessian(x, t, i = 0, j = 0, component = 0)
dx2_tt = dde.grad.hessian(x, t, i = 0, j = 0, component = 1)
pde1 = m1 * dx1_tt + k1 * (x_1 - L1) - k2 * (x_2 - x_1 - L1)
pde2 = m2 * dx2_tt + k2 * (x_2 - x_1 - L1)
return [pde1, pde2]
5.2.4. Define Initial ConditionsĀ¶
def boundary_init(t, on_boundary):
return on_boundary and np.isclose(t[0], 0)
geom = dde.geometry.Interval(0, t_max)
init_d1 = dde.icbc.PointSetBC(np.array([0]), np.array([d1_0]).reshape(-1, 1), component=0)
init_d2 = dde.icbc.PointSetBC(np.array([0]), np.array([d2_0]).reshape(-1, 1), component=1)
init_v1 = dde.OperatorBC(geom, lambda x, y, _: dy(x, y[:, [0]]), boundary_init)
init_v2 = dde.OperatorBC(geom, lambda x, y, _: dy(x, y[:, [1]]), boundary_init)
data = dde.data.PDE(geom,
pde,
[init_d1, init_d2, init_v1, init_v2],
num_domain = 2000,
num_boundary = 100,
num_test = 1000)
5.2.5. Define Network and Hyper-parametersĀ¶
layer_size = [1] + [20] * 3 + [2]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)
model = dde.Model(data, net)
model.compile("adam", lr = 1e-3)
5.2.6. Train & PredictionĀ¶
losshistory, train_state = model.train(epochs = 5000)
dde.saveplot(losshistory, train_state, issave = False, isplot = False)
dde.optimizers.config.set_LBFGS_options(maxiter=5000)
model.compile("L-BFGS")
losshistory, train_state = model.train()
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
5.2.7. EvaluationĀ¶
def vectorfield(w, t, p):
"""
Defines the differential equations for the coupled spring-mass system.
Arguments:
w : vector of the state variables:
w = [x1,y1,x2,y2]
t : time
p : vector of the parameters:
p = [m1,m2,k1,k2,L1,L2,b1,b2]
"""
x1, y1, x2, y2 = w
m1, m2, k1, k2, L1, L2, b1, b2 = p
# Create f = (x1',y1',x2',y2'):
f = [y1,
(-b1 * y1 - k1 * (x1 - L1) + k2 * (x2 - x1 - L2)) / m1,
y2,
(-b2 * y2 - k2 * (x2 - x1 - L2)) / m2]
return f
# Use ODEINT to solve the differential equations defined by the vector field
# Parameter values
# ODE solver parameters
abserr = 1.0e-8
relerr = 1.0e-6
stoptime = t_max
numpoints = 250
t = geom.random_points(5000)
t[:,0].sort()
# Create the time samples for the output of the ODE solver.
# I use a large number of points, only because I want to make
# a plot of the solution that looks nice.
# Pack up the parameters and initial conditions:
p = [m1, m2, k1, k2, L1, L2, b1, b2]
w0 = [d1_0, v1_0, d2_0, v2_0]
# Call the ODE solver.
wsol = odeint(vectorfield, w0, t[:, 0], args=(p,),
atol=abserr, rtol=relerr)
result = model.predict(t)
usol1 = np.array(result[:, 0])
usol2 = np.array(result[:, 1])
lw = 2
plt.figure(figsize = (12, 6))
plt.plot(t,
wsol[:, 0],
alpha = 1,
label = 'Mass1 GT',
c = 'r',
lw = lw
)
plt.plot(t,
wsol[:, 2],
alpha = 1,
label = 'Mass2 GT',
c = 'b',
lw = lw
)
plt.plot(t[:, 0],
usol1,
alpha = 1,
label = 'Mass1 Prediction',
linestyle='dashed',
c = 'k',
lw = lw
)
plt.plot(t[:, 0],
usol2,
alpha = 1,
label = 'Mass2 Prediction',
linestyle='dashed',
c = 'k',
lw = lw
)
plt.legend(fontsize = 15)
plt.xticks(fontsize = 25)
plt.yticks(fontsize = 25)
plt.xlabel('Time (s)', fontsize = 30)
plt.ylabel('Displacement (m)', fontsize = 30)
plt.show()
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')