2024.7.28周报

目录

摘要

ABSTRACT

一、文献阅读

一、题目

二、摘要

三、创新点

四、文献解读

一、Introduction

二、Saint-Venant方程

三、应用于水道建模的PINN

四、真实场景

五、结论

二、使用PINN解偏微分方程

摘要

本周阅读了一篇题目为Physics-Informed Neural Networks for Modeling Water Flows in a River Channel的论文,这篇论文提出了一种新的基于物理信息神经网络(PINN)的河道水流代理模型。本文研究PINN的性能直接从圣维南方程的配置建立在一个模拟的环境。这些实验揭示了水文建模的有希望的结果,并提出了替代方案,以解决传统方法中发现的主要挑战,同时帮助合成现实世界的表示。

ABSTRACT

This week, I read a paper titled "Physics-Informed Neural Networks for Modeling Water Flows in a River Channel." The paper proposes a new river channel water flow surrogate model based on Physics-Informed Neural Networks (PINNs). This study examines the performance of PINNs directly configured from the Saint-Venant equations in a simulated environment. These experiments reveal promising results for hydrological modeling and suggest alternatives to address the main challenges found in traditional methods while aiding in the synthesis of real-world representations.

一、文献阅读

一、题目

题目:Physics-Informed Neural Networks for Modeling Water Flows in a River Channel

链接:Physics-Informed Neural Networks for Modeling Water Flows in a River Channel | IEEE Journals & Magazine | IEEE Xplore

期刊\会议: ​​​​​​IEEE Transactions on Artificial Intelligence 

二、摘要

在实践中,防洪策略依赖于水文预报模型从概念或数据驱动的方法。受最近的研究成果的鼓舞,本研究提出了一种新的基于物理信息神经网络(PINNs)的河道水流代理模型。本文研究PINN的性能直接从圣维南方程的配置建立在一个模拟的环境。这些实验揭示了水文建模的有希望的结果,并提出了替代方案,以解决传统方法中发现的主要挑战,同时帮助合成现实世界的表示。在实际应用中,这些模型采用概念化或数据驱动的方法,概念模型要达到方法所提供的精度需要用到大量的物理参数.这些参数可能需要对环境有深入的了解,在复杂的盆地中可能很难识别。另一方面,虽然数据驱动的方法不需要这样的动态系统的知识,他们依赖于一个可靠的和有用的数据库,以保证系统行为的准确性。文章引入PINNs作为一个可行的解决方案,训练神经网络与一些训练数据和估计偏微分方程参数,管理潜在的动态。

In practice, flood control strategies rely on hydrological forecasting models from conceptual or data-driven approaches. Encouraged by the recent research results, this study proposes a new river flow agency model based on Physical Information Neural Networks (PINNs). This paper examines the performance of the PINN directly from the configuration of the San Venan program built in a simulated environment. These experiments reveal promising results for hydrological modeling and suggest alternatives to address major challenges found in traditional methods while aiding in the synthesis of real-world representations. In practice, these models adopt conceptual or data-driven methods, and a large number of physical parameters are required to achieve the accuracy provided by the methods. These parameters may require an in-depth understanding of the environment and can be difficult to identify in complex basins. On the other hand, while data-driven approaches do not require knowledge of such dynamic systems, they rely on a reliable and useful database to guarantee the accuracy of the system's behavior. This paper introduces PINNs as a feasible solution to train neural networks with some training data and estimate PDE parameters to manage potential dynamics.

三、创新点


(1)基于PINN的流域水流水文模拟新方法
(2)一种新的基于贝叶斯推理的PINN配置训练方法
(3)以圣维南方程为特征的综合区域水文模型分析

四、文献解读

一、Introduction

在文献中,发现了几种与洪水控制相关的水文预测策略,每种策略都以流域的特征和考虑的变量为特征。这些策略可以根据所采用的模型类型进行分类,可以是概念模型、经验模型或混合模型。基本上,概念模型是从所涉及的过程的物理推导出来的,因此,可以准确地描述潜在的现象。另一方面,经验模型采用启发式方法,并使用从观察中推断出的关系。这些数据驱动的模型不依赖于大量需要调整的物理参数,这些参数可能很难识别。最后,采用这两种策略的过程被称为混合建模。在此背景下,我们的工作提出了一种基于物理信息神经网络(pinn)的代理模型,用于建模通道流动。pinn是一种神经网络,它被训练来解决监督学习任务,同时尊重任何给定的物理定律。pinn的结构和训练策略使得吸收实际数据测量成为可能,同时应用物理正则化,从而导致控制潜在动力学的微分方程的解。除了微分方程的参数辨识之外,pinn在过程数据很少可用或随时间和空间分布不佳的情况下也很有吸引力。这种建模方法旨在规避概念和经验水文方法的缺点,同时保留其主要优点。

二、Saint-Venant方程

文献中用于建模明渠流的标准概念方法依赖于Saint-Venant方程的求解。两个偏微分方(PDEs)模型化了单一通道在没有下游流入情况下的一维动态。直观上,方程(1)和(2)是从质量守恒和动量守恒的物理定律中导出的。

三、应用于水道建模的PINN

一、让我们考虑由(1)和(2)给出的一维水流传播的动力学Q(x, t),边界和初始条件定义为:

在底面为2m、长度为2km、河床坡度为S_{0} = 0.0025、曼宁系数 n=0.02m^{-1/3}\cdot s的恒定矩形截面河道中。

二、物理信息神经网络的训练

为了训练神经网络NN_{Q}来模拟水流函数 Q(x, t) 并校准参数向量 λ = [S_{0}, n],需要一个适当的数据库来定义损失函数。训练数据U={<x_{u}^{i},t_{u}^{i},h_{u}^{i},Q_{u}^{i}>}_{i=1}^{N_{u}}可以通过给定的初始和边界条件从偏微分方程(PDE)的数值解得到。根据拉丁超立方抽样策略定义了一组配点C={<x_{c}^{i},t_{c}^{i},h_{c}^{i}>}_{i=1}^{N_{c}}

下图描述了对应于训练数据集的洪水波传播,其中N_{u}= 1200 个样本。请注意,训练包含了河道开始和结束时的河流动态,而正则化数据将包含内部点。这与典型的现实世界情景相似,其中河段的测量仅在战略点(即,感兴趣的城市)进行,这意味着运营商没有中间测量数据。

三、实验分析

为了研究PINN在模拟水流动态方面的性能和精确度,我们进行了实验,这些实验采用了不同数量的配点和不同的神经网络架构。(一)实验假设所有参数已知,专注于神经网络吸收测量数据的能力。(二)实验将数据吸收与系统特性化结合起来,假设Manning系数n和床面坡度S0是未知参数,需要识别。

损失函数为:

其中\psi =[x,t,h(x,t),\Delta h_{t}(x,t),\Delta h_{x}(x,t)],是由空间,时间,河流水位以及h(\cdot )的数值偏导数等变量组成的神经网络输入,输出是水流速率Q(x,t)。

为了最小化神经网络损失函数,使用500次ADAM算法,然后使用2500次有限记忆L-BFGS算法。ADAM算法是一种带动量的随机梯度下降算法,具有吸引人的收敛速度,此外,对于输入层,数据通过z-score方法进行标准化,神经网络的隐藏层由相同数量的N_{u}个神经元组成,每个神经元都具有双曲正切激活函数。表1显示了不同数量的隐藏层、每层神经元和搭配点的结果误差,而训练数据点的总数保持固定为Nu= 1200个样本。MSEu和MSEp分别是在搭配点和参数[S0, n] =[0.0025, 0.02]处P(NNQ, [S0, n]) = 0时的训练误差和由此产生的误差。为了量化神经网络近似流量函数Q(x, t)的能力,验证过程使用通道点x∈{500,1000,1500}m处的PDE的预测和解之间的均方误差指标。

下图展示了最佳情况下的这些结果,即表1中的实验2。图5(a)显示了损失函数作为训练历元的函数,而图5(b) - (d)分别给出了预测的时空解Q(x, t)与通道点x∈{500,1000,1500}m处的真实值的对
比。

四、真实场景

研究对象是巴西圣卡塔琳娜州大西洋沿岸的伊塔贾伊河流域。如下图所示,该流域是该州最大的,覆盖面积约15000平方公里。选定的实验区段是伊塔贾伊米林河,位于下图所示的矩形区域内,长约109公里。根据日本国际合作机构的研究,该地区的曼宁系数为0.032,河道的平均床面坡度为2.5996×10^{-4}。此外,该研究表明,分析区域的床面坡度变化在6.6667×10^{-5}5.8823×10^{-4}之间。Marcuzzo等人的研究也证实了这些发现,他们发现整个河道的平均床面坡度约为1.3×10^{-4}。实验的目的是证明所提出的方法能够找到这些参数,同时模型能够吸收真实数据。

该场景包括在河道入口和出口处的水位和流量测量。此外,政府机构提供了如横截面、面积和湿周的地理研究。假设通道条件规则,即路径没有显著变化,通过使用水位h和河道位置x的回归模型可以近似地估算出面积和湿周。在系统地分析了可用的地理信息后,面积和湿周的表达式由下式给出: 

注意,水位h是时间t和空间x的函数。由于函数A(x,h)和P_{w}(x,h)是线性的在参数\theta\alpha中,它们可以通过最小二乘算法有效地设置,满足在河道的初始、中间和终点的横截面动态。数据库用于训练PINN模型包括一个在每15分钟内测量的单一洪水事件,从2020年1月10日至14日,共有960个样本。数据图14描述了在入口点 x=0m和出口点 x=109000m的水位动态 h(x,t),以及洪水波的传播。下图报告了使用架构C和贝叶斯推理正则化的实验结果—公式。解决方案使用的初始均值向量 μ(0)=[0.00033,0.025]和方差 σ2(0)=[9×10−5,5×10−3],格式为 [S0,n]。下图(a) 显示了随着隐藏层的数量、每层的神经元数量和配点变化的MSE训练误差,而下图(b) 和 (c) 分别展示了每次实验中识别出的床面坡度 S0​ 和曼宁系数 n的参数。

五、结论

本文通过考虑PINN作为河道水流速率的替代水文模型,对目前的技术状况做出了贡献。该研究调查了PINN架构在直接从Saint-Venant方程配置实现的模拟环境中的性能。该配置包括一个受控的水文系统,具有恒定的横截面和河床坡度,以及没有测量噪声的合成数据。综合实验的结果令人鼓舞,表明pinn可以吸收数据并发现控制系统动力学的微分方程的参数。结果表明,PINN方法在不影响其主要优点的前提下,克服了概念模型参数估计和经验方法数据缺乏的缺点。因此,PINN方法是水文建模的可行替代方法。如果通过使用数据库的数值微分获得所需的导数,那么具有架构A的PINN与底层PDE模型的依附性最好。该模型预测水流的MSE验证误差在0.0263% ~ 0.305%之间,而曼宁系数达到n = 0.02176,河床斜率达到S0 = 0.002904,在最佳情况下接近真实值。另一方面,当沿河道的水文监测站数量有限时,架构C在建模过程中呈现出最有效的结果。为了创建架构C,我们提出了一种新的基于贝叶斯推理的pinn训练方法,当参数有较好的先验值时,得到了令
人满意的结果。

二、使用PINN解偏微分方程

网络架构:

import torch
import torch.nn as nn
import numpy as np
import matplotlib.pyplot as plt
from torch import autograd
import os
os.environ['KMP_DUPLICATE_LIB_OK']='True'
#使用最简单的全连接层
#用神经网络模拟微分方程,f'(x)=f(x),初始条件f(0)=1
class Net(nn.Module):
    #NL:n个l(线性,全连接)隐藏层,NN:输入数据的维数
    #NL就是有多少层隐藏层
    #NN是每层的神经元数量
    def __init__(self,NL,NN):
        super(Net,self).__init__()
        self.input_layer=nn.Linear(1,NN)
        self.hidden_layer=nn.Linear(NN,int(NN/2))
        self.output_layer=nn.Linear(int(NN/2),1)
    def forward(self,x):
        out=torch.tanh(self.input_layer(x))
        out=torch.tanh(self.hidden_layer(out))
        out_final=self.output_layer(out)
        return out_final
net=Net(4,20) #4层,20个
mse_cost_function=torch.nn.MSELoss(reduction='mean') #Mean squared error均方误差
optimizer=torch.optim.Adam(net.parameters(),lr=1e-4)#优化器

def ode_01(x,net):
    y=net(x)
    y_x=autograd.grad(y,x,grad_outputs=torch.ones_like(net(x)),create_graph=True)[0]
    return y-y_x   #y-y'=0
plt.ion()
iterations=200000
for epoch in range(iterations):
    optimizer.zero_grad()
    #求边界条件的损失函数
    x_0=torch.zeros(2000,1)
    y_0=net(x_0)
    mse_i=mse_cost_function(y_0,torch.ones(2000,1)) #f(0)-1=0
    #方程的损失函数
    x_in=np.random.uniform(low=0.0,high=2.0,size=(2000,1))
    pt_x_in=autograd.Variable(torch.from_numpy(x_in).float(),requires_grad=True)
    pt_y_colection=ode_01(pt_x_in,net)
    pt_all_zeros=autograd.Variable(torch.from_numpy(np.zeros((2000,1))).float(),requires_grad=False)
    mse_f=mse_cost_function(pt_y_colection,pt_all_zeros) #y-y'=0

    loss=mse_i+mse_f
    loss.backward()
    optimizer.step()

    if epoch%1000==0:
        y=torch.exp(pt_x_in) #y 真实值
        y_train0=net(pt_x_in)
        print(epoch,"Training Loss",loss.data)
        print(f'times{epoch} - loss:{loss.item()} -y_0:{y_0}')
        plt.cla()
        plt.scatter(pt_x_in.detach().numpy(),y.detach().numpy())
        plt.scatter(pt_x_in.detach().numpy(),y_train0.detach().numpy(),c='red')
        plt.pause(0.1)

实验:

"""
A scratch for PINN solving the following PDE
u_xx-u_yyyy=(2-x^2)*exp(-y)
Author: ST
Date: 2023/2/26
"""
import torch
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import os
os.environ['KMP_DUPLICATE_LIB_OK']='True'
epochs = 1000    # 训练代数
h = 100    # 画图网格密度
N = 1000    # 内点配置点数
N1 = 100    # 边界点配置点数
N2 = 1000    # PDE数据点

def setup_seed(seed):
    torch.manual_seed(seed)
    torch.cuda.manual_seed_all(seed)
    torch.backends.cudnn.deterministic = True

# 设置随机数种子
setup_seed(888888)

# Domain and Sampling
def interior(n=N):
    # 内点
    x = torch.rand(n, 1)
    y = torch.rand(n, 1)
    cond = (2 - x ** 2) * torch.exp(-y)
    return x.requires_grad_(True), y.requires_grad_(True), cond


def down_yy(n=N1):
    # 边界 u_yy(x,0)=x^2
    x = torch.rand(n, 1)
    y = torch.zeros_like(x)
    cond = x ** 2
    return x.requires_grad_(True), y.requires_grad_(True), cond


def up_yy(n=N1):
    # 边界 u_yy(x,1)=x^2/e
    x = torch.rand(n, 1)
    y = torch.ones_like(x)
    cond = x ** 2 / torch.e
    return x.requires_grad_(True), y.requires_grad_(True), cond


def down(n=N1):
    # 边界 u(x,0)=x^2
    x = torch.rand(n, 1)
    y = torch.zeros_like(x)
    cond = x ** 2
    return x.requires_grad_(True), y.requires_grad_(True), cond


def up(n=N1):
    # 边界 u(x,1)=x^2/e
    x = torch.rand(n, 1)
    y = torch.ones_like(x)
    cond = x ** 2 / torch.e
    return x.requires_grad_(True), y.requires_grad_(True), cond


def left(n=N1):
    # 边界 u(0,y)=0
    y = torch.rand(n, 1)
    x = torch.zeros_like(y)
    cond = torch.zeros_like(x)
    return x.requires_grad_(True), y.requires_grad_(True), cond


def right(n=N1):
    # 边界 u(1,y)=e^(-y)
    y = torch.rand(n, 1)
    x = torch.ones_like(y)
    cond = torch.exp(-y)
    return x.requires_grad_(True), y.requires_grad_(True), cond

def data_interior(n=N2):
    # 内点
    x = torch.rand(n, 1)
    y = torch.rand(n, 1)
    cond = (x ** 2) * torch.exp(-y)
    return x.requires_grad_(True), y.requires_grad_(True), cond


# Neural Network
class MLP(torch.nn.Module):
    def __init__(self):
        super(MLP, self).__init__()
        self.net = torch.nn.Sequential(
            torch.nn.Linear(2, 32),
            torch.nn.Tanh(),
            torch.nn.Linear(32, 32),
            torch.nn.Tanh(),
            torch.nn.Linear(32, 32),
            torch.nn.Tanh(),
            torch.nn.Linear(32, 32),
            torch.nn.Tanh(),
            torch.nn.Linear(32, 1)
        )

    def forward(self, x):
        return self.net(x)


# Loss
loss = torch.nn.MSELoss()


def gradients(u, x, order=1):
    if order == 1:
        return torch.autograd.grad(u, x, grad_outputs=torch.ones_like(u),
                                   create_graph=True,
                                   only_inputs=True, )[0]
    else:
        return gradients(gradients(u, x), x, order=order - 1)

# 以下7个损失是PDE损失
def l_interior(u):
    # 损失函数L1
    x, y, cond = interior()
    uxy = u(torch.cat([x, y], dim=1))
    return loss(gradients(uxy, x, 2) - gradients(uxy, y, 4), cond)


def l_down_yy(u):
    # 损失函数L2
    x, y, cond = down_yy()
    uxy = u(torch.cat([x, y], dim=1))
    return loss(gradients(uxy, y, 2), cond)


def l_up_yy(u):
    # 损失函数L3
    x, y, cond = up_yy()
    uxy = u(torch.cat([x, y], dim=1))
    return loss(gradients(uxy, y, 2), cond)


def l_down(u):
    # 损失函数L4
    x, y, cond = down()
    uxy = u(torch.cat([x, y], dim=1))
    return loss(uxy, cond)


def l_up(u):
    # 损失函数L5
    x, y, cond = up()
    uxy = u(torch.cat([x, y], dim=1))
    return loss(uxy, cond)


def l_left(u):
    # 损失函数L6
    x, y, cond = left()
    uxy = u(torch.cat([x, y], dim=1))
    return loss(uxy, cond)


def l_right(u):
    # 损失函数L7
    x, y, cond = right()
    uxy = u(torch.cat([x, y], dim=1))
    return loss(uxy, cond)

# 构造数据损失
def l_data(u):
    # 损失函数L8
    x, y, cond = data_interior()
    uxy = u(torch.cat([x, y], dim=1))
    return loss(uxy, cond)


# Training

u = MLP()
opt = torch.optim.Adam(params=u.parameters())

for i in range(epochs):
    opt.zero_grad()
    l = l_interior(u) \
        + l_up_yy(u) \
        + l_down_yy(u) \
        + l_up(u) \
        + l_down(u) \
        + l_left(u) \
        + l_right(u) \
        + l_data(u)
    l.backward()
    opt.step()
    if i % 100 == 0:
        print(i)

# Inference
xc = torch.linspace(0, 1, h)
xm, ym = torch.meshgrid(xc, xc,indexing='ij')
xx = xm.reshape(-1, 1)
yy = ym.reshape(-1, 1)
xy = torch.cat([xx, yy], dim=1)
u_pred = u(xy)
u_real = xx * xx * torch.exp(-yy)
u_error = torch.abs(u_pred-u_real)
u_pred_fig = u_pred.reshape(h,h)
u_real_fig = u_real.reshape(h,h)
u_error_fig = u_error.reshape(h,h)
print("Max abs error is: ", float(torch.max(torch.abs(u_pred - xx * xx * torch.exp(-yy)))))
# 仅有PDE损失    Max abs error:  0.004852950572967529
# 带有数据点损失  Max abs error:  0.0018916130065917969

# 作PINN数值解图
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(xm.detach().numpy(), ym.detach().numpy(), u_pred_fig.detach().numpy())
ax.text2D(0.5, 0.9, "PINN", transform=ax.transAxes)
plt.show()
fig.savefig("PINN solve.png")

# 作真解图
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(xm.detach().numpy(), ym.detach().numpy(), u_real_fig.detach().numpy())
ax.text2D(0.5, 0.9, "real solve", transform=ax.transAxes)
plt.show()
fig.savefig("real solve.png")

# 误差图
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(xm.detach().numpy(), ym.detach().numpy(), u_error_fig.detach().numpy())
ax.text2D(0.5, 0.9, "abs error", transform=ax.transAxes)
plt.show()
plt.pause(0.1)
fig.savefig("abs error.png")
print("u_pred_fig shape:", u_pred_fig.shape)
print("Data range:", u_pred_fig.min(), u_pred_fig.max())

 

实验结果:

  • 18
    点赞
  • 28
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值