2025.2.23机器学习笔记:PINN文献阅读

news2025/2/25 13:55:17

2025.2.23周报

  • 一、文献阅读
    • 题目信息
    • 摘要
    • Abstract
    • 创新点
    • 网络架构
      • 架构A
      • 架构B
      • 架构C
    • 实验
    • 结论
    • 后续展望

一、文献阅读

题目信息

  • 题目: Physics-Informed Neural Networks for Modeling Water Flows in a River Channel
  • 期刊: IEEE TRANSACTIONS ON ARTIFICIAL INTELLIGENCE
  • 作者: Luis Fernando Nazari,Eduardo Camponogara,Laio Oriel Seman
  • 发表时间: 2024/03/03
  • 文章链接: https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=9863669

摘要

洪水灾害对世界各地的影响重大,造成了严重的社会和经济问题,洪水防控策略依赖于水文预报模型,而这些模型主要通过数据驱动方法获得。所以模型需要达到期望的精度需要大量物理参数,而这些参数往往难以确获取;且传统的数据驱动模型往往需要可靠且真实的数据保证模型的数值模拟准确性。所以在以上背景研究下,本论文提出了一个基于物理信息神经网络(Physics-Informed Neural Networks,PINNs)的河道水流替代模型,作者总共设计了三个PINNs模型架构A、B、C,通过调节模型的输入与输出参数以探究在不同场景中模型的适用性。其中在架构C中作者引入了贝叶斯推断,限制参数偏离合理范围,提升物理可解释性。该论文旨在解决传统方法存在的问题,为流域水流的水文建模领域的研究做出贡献。

Abstract

Flood disasters exert profound impacts on populations worldwide, resulting in severe social and economic challenges. Effective flood prevention strategies hinge on hydrological forecasting models, which are predominantly derived through data-driven approaches. Achieving the desired precision in these models necessitates a multitude of physical parameters, which are often difficult to accurately acquire. Moreover, conventional data-driven models typically require reliable and authentic datasets to ensure the accuracy of numerical simulations. Against this backdrop, this paper proposes a surrogate model for river channel water flow based on Physics-Informed Neural Networks (PINNs). The authors developed three distinct PINN architectures—A, B, and C—by adjusting the input and output parameters to investigate their applicability across varying scenarios. Notably, in Architecture C, Bayesian inference is introduced to constrain parameter deviations within physically plausible ranges, thereby enhancing physical interpretability. This study aims to address the limitations inherent in traditional methodologies, contributing to advancements in the field of hydrological modeling for watershed water flows.

创新点

  1. 将PINNs应用于河道水流建模,通过融合圣维南方程的物理约束与实测数据,解决了传统方法中数据稀缺和参数依赖的问题;
  2. 通过PINNs直接预测水文动力学的关键参数,如:曼宁系数 n n n和河床坡度 S 0 S_0 S0。这种方法减少了传统方法中对地区数据的依赖;
  3. 使用了基于贝叶斯推理的PINNs训练新方法,解决参数预测的不确定性问题。

网络架构

本文提出了一种基于物理信息神经网络(PINN)的河流水流建模方法,重点解决了传统水文模型中参数识别困难和数据不足的问题。
作者提出了圣维南方程作为物理模型:
在这里插入图片描述
关于 Q ( x , t ) Q(x,t) Q(x,t) P w ( x , t ) P_w(x,t) Pw(x,t)的形象表示如下图所示:
在这里插入图片描述
下面对论文中三种网络架构进行分析:

架构A

  1. 输入
    ψ = [ x , t , h ( x , t ) , Δ h t ( x , t ) , Δ h x ( x , t ) ] ∈ R 5 \psi = [x, t, h(x,t), \Delta h_t(x,t), \Delta h_x(x,t)] \in \mathbb{R}^5 ψ=[x,t,h(x,t),Δht(x,t),Δhx(x,t)]R5
    其中,参数含义如下:
    x x x:空间;
    t t t:时间;
    h ( x , t ) h(x,t) h(x,t):水位;
    Δ h t \Delta h_t Δht
    数值微分计算的时间导数;
    Δ h x \Delta h_x Δhx:空间导数。
  2. 输出
    模型的输出为 Q Q Q水流速。
  3. 损失函数
    PINN的模型架构的损失函数由两部分组成:数据误差物理误差,即 Loss = MSE u + MSE p \text{Loss} = \text{MSE}_u + \text{MSE}_p Loss=MSEu+MSEp
    M S E u MSE_u MSEu为数据误差,公式表达为 MSE u = 1 N u ∑ i = 1 N u ∣ N N Q ( ψ u i ) − Q u i ∣ 2 \text{MSE}_u = \frac{1}{N_u} \sum_{i=1}^{N_u} \left| \mathcal{NN}_Q(\psi_u^i) - Q_u^i \right|^2 MSEu=Nu1i=1Nu NNQ(ψui)Qui 2,表示真实值与预测值之间的误差。具体展开为: MSE p = 1 N c ∑ i = 1 N c [ ∣ ∂ A ∂ t + ∂ N N Q ∂ x ∣ 2 + ∣ ∂ N N Q ∂ t + ∂ ∂ x ( N N Q 2 A ) + g A ( ∂ h ∂ x + S f − S 0 ) ∣ 2 ] \text{MSE}_p = \frac{1}{N_c} \sum_{i=1}^{N_c} \left[ \left| \frac{\partial A}{\partial t} + \frac{\partial \mathcal{NN}_Q}{\partial x} \right|^2 + \left| \frac{\partial \mathcal{NN}_Q}{\partial t} + \frac{\partial}{\partial x} \left( \frac{\mathcal{NN}_Q^2}{A} \right) + gA \left( \frac{\partial h}{\partial x} + S_f - S_0 \right) \right|^2 \right] MSEp=Nc1i=1Nc[ tA+xNNQ 2+ tNNQ+x(ANNQ2)+gA(xh+SfS0) 2]

在这里插入图片描述

架构B

  1. input
    ψ = [ x , t , h ( x , t ) ] ∈ R 3 \psi = [x, t, h(x,t)] \in \mathbb{R}^3 ψ=[x,t,h(x,t)]R3
    其中,参数含义如下:
    x x x:空间;
    t t t:时间;
    h ( x , t ) h(x,t) h(x,t):水位程高;
  2. output
    输出为: h和Q,分别表示水位高度与流速。
  3. Loss function
    损失函数扩展为三部分,分别为:数据误差圣维南方程质量守恒约束以及圣维南方程动量守恒约束
    l o s s f u n c t i o n = 1 N u ∑ i = 1 N u ∣ N N ( ψ u i ) − ( h i , Q i ) ∣ 2 + 1 N c ∑ i = 1 N c ∣ ∂ A ∂ h ∂ ∂ t N N h ( ψ c i ) + ∂ ∂ x N N Q ( ψ c i ) ∣ 2 + 1 N c ∑ i = 1 N c ∣   ∂ ∂ t N N Q ( ψ c i ) + ∂ ∂ x N N Q 2 ( ψ c i ) A ( x c i , t c i ) + g A ( x c i , t c i ) ( ∂ ∂ x N N h ( ψ c i ) + S f ( x c i , t c i ) − S 0 ) ∣ 2 loss function = \frac{1}{N_{u}} \sum_{i=1}^{N_{u}}\left|\mathcal{N N}\left(\psi_{u}^{i}\right)-\left(h^{i}, Q^{i}\right)\right|^{2} +\frac{1}{N_{c}} \sum_{i=1}^{N_{c}}\left|\frac{\partial A}{\partial h} \frac{\partial}{\partial t} \mathcal{N} \mathcal{N}_{h}\left(\psi_{c}^{i}\right)+\frac{\partial}{\partial x} \mathcal{N} \mathcal{N}_{Q}\left(\psi_{c}^{i}\right)\right|^{2} +\frac{1}{N_{c}} \sum_{i=1}^{N_{c}} \left\lvert\, \frac{\partial}{\partial t} \mathcal{N N}_{Q}\left(\psi_{c}^{i}\right)+\frac{\partial}{\partial x} \frac{\mathcal{N N}_{Q}^{2}\left(\psi_{c}^{i}\right)}{A\left(x_{c}^{i}, t_{c}^{i}\right)}\right. +\left.g A\left(x_{c}^{i}, t_{c}^{i}\right)\left(\frac{\partial}{\partial x} \mathcal{N} \mathcal{N}_{h}\left(\psi_{c}^{i}\right)+S_{f}\left(x_{c}^{i}, t_{c}^{i}\right)-S_{0}\right)\right|^{2} lossfunction=Nu1i=1Nu NN(ψui)(hi,Qi) 2+Nc1i=1Nc hAtNNh(ψci)+xNNQ(ψci) 2+Nc1i=1Nc tNNQ(ψci)+xA(xci,tci)NNQ2(ψci)+gA(xci,tci)(xNNh(ψci)+Sf(xci,tci)S0) 2
    其中, N c N_c Nc为配准点(collocation points)的数量,这些点是用来强制满足物理方程的时空坐标; N N h ( ψ c i ) \mathcal{N} \mathcal{N}_{h}\left(\psi_{c}^{i}\right) NNh(ψci):神经网络预测的水位高度h在配准点 ( x c i , t c i ) \left(x_c^i, t_c^i\right) (xci,tci)上的值; N N Q ( ψ c i ) \mathcal{N} \mathcal{N}_{Q}\left(\psi_{c}^{i}\right) NNQ(ψci):神经网络预测的水流量Q在配准 ( x c i , t c i ) \left(x_c^i, t_c^i\right) (xci,tci)上的值。
    在这里插入图片描述

架构C

  1. input
    ψ = [ x , t , h ( x , t ) , Δ h t ( x , t ) ] ∈ R 4 \psi = [x, t, h(x,t), \Delta h_t(x,t)] \in \mathbb{R}^4 ψ=[x,t,h(x,t),Δht(x,t)]R4
    其中,参数 R 4 \mathbb{R}^4 R4表示四维空间;
    x x x:空间;
    t t t:时间;
    h h h:水位;
    Δ h t \Delta h_t Δht:表示水面程高对时间t的求导。
  2. 输出
    输出与架构B相同都是 h h h Q Q Q,输出是一个二维向量。
  3. 损失函数
    与架构A和架构B不同的是,架构C在损失函数引入贝叶斯推断,增加惩罚项,其整体的损失函数表示为:
    Loss Function = 1 N u ∑ i = 1 N u ∣ N N ( ψ u i ) − [ h i , Q i ] ∣ 2 ⏟ 数据误差 + 1 N c ∑ i = 1 N c ∣ ∂ A ∂ h ∂ N N h ∂ t + ∂ N N Q ∂ x ∣ 2 ⏟ 质量守恒 + 1 N c ∑ i = 1 N c ∣ ∂ N N Q ∂ t + ∂ ∂ x N N Q 2 A + g A ( ∂ N N h ∂ x + S f − S 0 ) ∣ 2 ⏟ 动量守恒 + ∑ i = 1 N p 1 σ i 2 ( λ i − μ i ) 2 ⏟ 贝叶斯正则化 + 1 N c ∑ i = 1 N c ∣ ∂ N N h ∂ t − Δ h t ∣ 2 ⏟ 导数误差 \text{Loss Function} = \underbrace{\frac{1}{N_u} \sum_{i=1}^{N_u} \left| \mathcal{NN}(\psi_u^i) - [h^i, Q^i] \right|^2}_{\text{数据误差}} + \underbrace{\frac{1}{N_c} \sum_{i=1}^{N_c} \left| \frac{\partial A}{\partial h} \frac{\partial \mathcal{NN}_h}{\partial t} + \frac{\partial \mathcal{NN}_Q}{\partial x} \right|^2}_{\text{质量守恒}} + \underbrace{\frac{1}{N_c} \sum_{i=1}^{N_c} \left| \frac{\partial \mathcal{NN}_Q}{\partial t} + \frac{\partial}{\partial x} \frac{\mathcal{NN}_Q^2}{A} + g A \left( \frac{\partial \mathcal{NN}_h}{\partial x} + S_f - S_0 \right) \right|^2}_{\text{动量守恒}} + \underbrace{\sum_{i=1}^{N_p} \frac{1}{\sigma_i^2} \left( \lambda_i - \mu_i \right)^2}_{\text{贝叶斯正则化}} + \underbrace{\frac{1}{N_c} \sum_{i=1}^{N_c} \left| \frac{\partial \mathcal{NN}_h}{\partial t} - \Delta h_t \right|^2}_{\text{导数误差}} Loss Function=数据误差 Nu1i=1Nu NN(ψui)[hi,Qi] 2+质量守恒 Nc1i=1Nc hAtNNh+xNNQ 2+动量守恒 Nc1i=1Nc tNNQ+xANNQ2+gA(xNNh+SfS0) 2+贝叶斯正则化 i=1Npσi21(λiμi)2+导数误差 Nc1i=1Nc tNNhΔht 2
    引入贝叶斯推断可以参数 S 0 S_0 S0 n n n 更靠,贝叶斯推断表示为: 1 N p ∑ j = 1 N p 1 σ i 2 ( λ i − μ i ) 2 \frac{1}{N_p} \sum_{j=1}^{N_p} \frac{1}{\sigma_i^2} (\lambda_i - \mu_i)^2 Np1j=1Npσi21(λiμi)2
    其中, λ i \lambda_i λi为待估计参数; N p N_p Np为参数数量; μ i \mu_i μi参数的均值; σ i 2 {\sigma_i^2} σi2为参数的方差。
    在这里插入图片描述

架构A、B、C的比较以及适用场景

架构对输入参数的依赖车程度物理可解释性适用场景
A密集监测场景
B监测站稀疏的实际场景
C最高需平衡数据与物理定律的场景

实验

  1. PINNs用于河道水流建模的初步实验结果
    在模拟环境下,以圣维南方程配置构建的系统中进行实验,目的是在评估PINNs预测河道水流的有效性。实验分两个阶段,第一阶段假设PDE参数和水位函数已知,验证PINNs的数据同化能力;第二阶段聚焦于识别定义PDE的向量参数,以描述水流动力学并捕捉河道内在特征。
    通过数值解圣维南方程得到训练数据和配置点集,采用如架构A进行实验,结果表明PINN方法在系统识别方面,使用相对少的数据点,在所有实验中达到相当的精度和强大的数据同化能力,验证的均方误差在0.00648% - 0.0466%之间。
    在这里插入图片描述

  2. 不同架构下的建模结果
    架构A
    在将 S 0 S_0 S0 n n n作为优化变量加入损失函数后,经过一定训练轮次,架构A在所有实验中的预测误差(MSE)低于0.305%,最佳情况MSE为0.0263%。
    在这里插入图片描述
    基于PINNs发现的摩擦斜率和曼宁系数近似于圣维南方程中的实际参数,表明PINNs通过利用微分方程模型知识实现有效正则化。尽管没有理论保证神经网络模型和参数发现的联合学习收敛到全局最优,但在一定条件下可达到较好预测精度。
    在这里插入图片描述
    架构B
    为解决数据点缺乏问题提出架构B,其输入变量仅包含距离、时间和河流水位,输出除水流率外还尝试学习水位函数用于自动微分计算物理正则化中的偏导数。
    在这里插入图片描述
    实验结果显示架构B发现真实PDE参数的能力不如架构A,但能较好地近似系统行为,MSE验证指数较低。
    在这里插入图片描述
    架构C
    为改进参数识别提出架构C,采用新结构并受益于新的损失函数,引入贝叶斯推理到训练过程。
    架构C添加 Δ h t Δh_t Δht作为输入变量,新的损失函数包含对参数偏离均值的惩罚项和对输出自动微分与数值导数差异的惩罚项。
    实验结果表明,架构C在系统行为的泛化方面,MSE验证指数在0.0762% - 0.298%之间,相比架构B有所改进,在参数估计方面也有显著提升,能以较好的初始先验值逼近真实参数并保持物理特性。对比固定方差向量的情况,虽然参数有显著改进,但会导致系统识别能力的损失。
    在这里插入图片描述

  3. 伊塔雅伊米林河中的实验结果
    以伊塔雅伊河的伊塔雅伊米林河河段为研究对象,该区域有水位和流量测量数据,以及政府提供的地理研究数据。
    在这里插入图片描述
    利用2020年1月10 - 14日一次洪水事件的960个样本数据训练PINN模型,采用架构C并加入贝叶斯推理正则化的实验结果显示,架构为八层隐藏层、每层30个神经元、80000个配置点时,数据集同化效果最佳(MSE = 4.126×10⁻³),发现的床坡和曼宁系数参数,其中曼宁系数收敛到接近文献中的值,床坡参数收敛到区域研究中的极限值内。
    在这里插入图片描述
    在这里插入图片描述
    代码如下:

# 导入必要的库
import torch  # 用于张量计算和神经网络搭建
import torch.nn as nn  # 提供神经网络层和激活函数
import numpy as np  # 用于数据转换和绘图
import matplotlib.pyplot as plt  # 用于可视化训练结果
from time import time  # 用于测量训练时间
from pytorch_lbfgs import LBFGS  # 导入支持GPU的L-BFGS优化器

# 设置随机种子,确保结果可重复
torch.manual_seed(42)

# 检查并设置计算设备,确保使用GPU
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
if not torch.cuda.is_available():
    raise RuntimeError("GPU not available. This code requires CUDA support.")
print(f"Using device: {device} ({torch.cuda.get_device_name(0)})")


# 定义Saint-Venant方程
def saint_venant_equations(h, Q, x, t, S0, n, g=9.81, eps=1e-6):
    """计算Saint-Venant方程的残差,用于物理约束"""
    width = 2.0  # 矩形横截面宽度2m
    A = h * width + eps  # 横截面积,添加偏移避免除零
    P_w = 2 * h + width + eps  # 湿周长,添加偏移
    R = A / P_w  # 水力半径

    # 自动微分计算偏导数
    h_t = torch.autograd.grad(h.sum(), t, create_graph=True)  # 计算 ∂h/∂t
    h_t = h_t[0] if h_t is not None else torch.zeros_like(h)  # 若梯度为空,返回零张量
    h_x = torch.autograd.grad(h.sum(), x, create_graph=True)  # 计算 ∂h/∂x
    h_x = h_x[0] if h_x is not None else torch.zeros_like(h)  # 若梯度为空,返回零张量
    Q_x = torch.autograd.grad(Q.sum(), x, create_graph=True)  # 计算 ∂Q/∂x
    Q_x = Q_x[0] if Q_x is not None else torch.zeros_like(Q)  # 若梯度为空,返回零张量
    Q_t = torch.autograd.grad(Q.sum(), t, create_graph=True)  # 计算 ∂Q/∂t
    Q_t = Q_t[0] if Q_t is not None else torch.zeros_like(Q)  # 若梯度为空,返回零张量
    Q2_A = Q ** 2 / (A + eps)  # 计算 Q^2/A,添加偏移避免除零
    Q2_A_x = torch.autograd.grad(Q2_A.sum(), x, create_graph=True)  # 计算 ∂(Q^2/A)/∂x
    Q2_A_x = Q2_A_x[0] if Q2_A_x is not None else torch.zeros_like(Q)  # 若梯度为空,返回零张量

    S_f = (n ** 2 * Q * torch.abs(Q)) / (A ** 2 * R ** (4 / 3) + eps)  # 摩擦坡度,添加偏移
    f1 = (width * h_t) + Q_x  # 质量守恒方程
    f2 = Q_t + Q2_A_x + g * A * (h_x + S_f - S0)  # 动量守恒方程

    return f1, f2  # 返回残差


# 数据生成函数
def generate_data(N_u=1200, N_c=150000, N_val=100):
    """生成训练数据、配准点和验证点,模拟2km河道洪水波"""
    t = torch.linspace(0, 12000, N_u // 2, device=device, dtype=torch.float32).reshape(-1, 1)  # 生成时间序列
    x_train = torch.cat([torch.zeros_like(t), 2000 * torch.ones_like(t)], dim=0)  # 设置x=0m和x=2000m
    t_train = torch.cat([t, t], dim=0)  # 拼接完整时间序列

    h_train = torch.zeros_like(t_train)  # 初始化水深数组
    mask1 = t_train <= 6000  # 上升阶段
    mask2 = (t_train > 6000) & (t_train <= 9000)  # 下降阶段
    mask3 = t_train > 9000  # 平稳阶段
    h_train[mask1] = (1.000735) ** (t_train[mask1] / 10) * 0.4  # 上升阶段水深
    h_ta = (1.000735) ** (6000 / 10) * 0.4  # t=6000s参考值
    h_train[mask2] = (0.9985) ** ((t_train[mask2] - 6000) / 10) * h_ta  # 下降阶段水深
    h_train[mask3] = 0.4  # 平稳阶段水深

    A0 = h_train * 2.0 + 1e-6  # 初始横截面积
    P0 = 2 * h_train + 2.0 + 1e-6  # 初始湿周长
    Q_train = ((0.0025) ** 0.5 / 0.02) * (A0 ** (5 / 3)) / (P0 ** (2 / 3))  # 稳态流量

    x_c = torch.rand(N_c, 1, device=device, dtype=torch.float32) * 2000  # 配准点空间
    t_c = torch.rand(N_c, 1, device=device, dtype=torch.float32) * 12000  # 配准点时间

    x_val = torch.tensor([500, 1000, 1500], device=device, dtype=torch.float32).repeat(N_val).reshape(-1, 1)  # 验证点位置
    t_val = torch.linspace(0, 12000, N_val, device=device, dtype=torch.float32).repeat(3).reshape(-1, 1)  # 验证点时间

    return (x_train.requires_grad_(True), t_train.requires_grad_(True), h_train, Q_train), \
        (x_c.requires_grad_(True), t_c.requires_grad_(True)), \
        (x_val.requires_grad_(True), t_val.requires_grad_(True))


# 定义PINN模型(架构C)
class PINN_C(nn.Module):
    def __init__(self, hidden_layers=8, neurons=20):
        """定义架构C的神经网络:输入[x, t, Δh_t],输出[h, Q]"""
        super(PINN_C, self).__init__()
        layers = [nn.Linear(3, neurons), nn.Tanh()]  # 输入层:3维 -> neurons,使用tanh激活
        for _ in range(hidden_layers - 1):
            layers += [nn.Linear(neurons, neurons), nn.Tanh()]  # 添加隐层
        layers += [nn.Linear(neurons, 2)]  # 输出层:neurons -> 2 (h, Q)
        self.net = nn.Sequential(*layers).to(device)  # 构建网络并移到GPU
        self.S0 = nn.Parameter(torch.tensor([0.002], device=device, dtype=torch.float32))  # 初始化床坡
        self.n = nn.Parameter(torch.tensor([0.03], device=device, dtype=torch.float32))  # 初始化曼宁系数

    def forward(self, x, t, dh_dt):
        """前向传播,输入[x, t, Δh_t],输出[h, Q]"""
        inputs = torch.cat([x, t, dh_dt], dim=1)  # 拼接输入向量
        out = self.net(inputs)  # 通过神经网络
        h = torch.relu(out[:, 0:1]) + 0.01  # 输出水深,确保非负并加偏移
        Q = torch.relu(out[:, 1:2]) + 0.01  # 输出流量,确保非负并加偏移
        return h, Q


# 损失函数
def loss_fn_C(model, x_u, t_u, h_u, Q_u, x_c, t_c, mu, sigma):
    """计算损失函数,包括数据项、物理项、正则化项和导数惩罚项"""
    dh_dt_u = torch.zeros_like(h_u)  # 初始化训练数据的 Δh_t
    dh_dt_u[1:] = (h_u[1:] - h_u[:-1]) / 10  # 计算数值导数
    dh_dt_u[0] = dh_dt_u[1]  # 边界填充

    h_pred, Q_pred = model(x_u, t_u, dh_dt_u)  # 预测训练点的 h 和 Q
    mse_u = torch.mean((h_pred - h_u) ** 2 + (Q_pred - Q_u) ** 2)  # 计算数据误差

    dh_dt_c = torch.zeros_like(x_c)  # 初始化配准点的 Δh_t
    h_c, Q_c = model(x_c, t_c, dh_dt_c)  # 预测配准点的 h 和 Q
    f1, f2 = saint_venant_equations(h_c, Q_c, x_c, t_c, model.S0, model.n)  # 计算PDE残差
    mse_p = torch.mean(f1 ** 2) + torch.mean(f2 ** 2)  # 计算物理误差

    reg = (1 / (sigma[0] ** 2 + 1e-6)) * (model.S0 - mu[0]) ** 2 + \
          (1 / (sigma[1] ** 2 + 1e-6)) * (model.n - mu[1]) ** 2  # 参数正则化项

    h_t = torch.autograd.grad(h_c.sum(), t_c, create_graph=True)  # 计算配准点 ∂h/∂t
    h_t = h_t[0] if h_t is not None else torch.zeros_like(h_c)  # 若梯度为空,返回零张量
    penalty = torch.mean((h_t - dh_dt_c) ** 2)  # 计算导数惩罚项

    return mse_u + mse_p + reg + penalty  # 返回总损失


# 验证MSE计算
def compute_validation_mse(model, x_val, t_val, h_true, Q_true):
    """计算验证点的MSE"""
    dh_dt_val = torch.zeros_like(h_true)  # 初始化验证点 Δh_t
    dh_dt_val[1:] = (h_true[1:] - h_true[:-1]) / 10  # 计算数值导数
    dh_dt_val[0] = dh_dt_val[1]  # 边界填充
    with torch.no_grad():  # 无梯度计算
        h_pred, Q_pred = model(x_val, t_val, dh_dt_val)  # 预测验证点 h 和 Q
        mse_val = torch.mean((h_pred - h_true) ** 2 + (Q_pred - Q_true) ** 2)  # 计算MSE
    return mse_val.item()  # 返回标量值


# 训练函数
def train_model(model, x_u, t_u, h_u, Q_u, x_c, t_c, x_val, t_val, h_val, Q_val):
    """训练模型:2000次ADAM + 5000次L-BFGS,全部在GPU上"""
    optimizer_adam = torch.optim.Adam(model.parameters(), lr=0.0001)  # 初始化ADAM优化器
    optimizer_lbfgs = LBFGS(model.parameters(), lr=0.1, max_iter=20)  # 初始化GPU支持的L-BFGS

    mu = torch.tensor([0.002, 0.03], device=device, dtype=torch.float32)  # 初始化先验均值
    sigma = torch.tensor([4.17e-4, 5e-3], device=device, dtype=torch.float32)  # 初始化先验方差

    S0_history, n_history, losses, val_mses = [], [], [], []  # 初始化历史记录列表

    print("Training with ADAM on GPU (2000 epochs)...")
    start_time = time()  # 记录ADAM阶段开始时间
    for epoch in range(2000):  # 2000次ADAM迭代
        optimizer_adam.zero_grad()  # 清零梯度
        loss = loss_fn_C(model, x_u, t_u, h_u, Q_u, x_c, t_c, mu, sigma)  # 计算总损失
        if torch.isnan(loss):  # 检查是否出现NaN
            print(f"NaN detected at epoch {epoch}")
            break
        loss.backward()  # 反向传播
        torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm=1.0)  # 梯度裁剪
        optimizer_adam.step()  # 更新参数

        S0_history.append(model.S0.item())  # 记录床坡
        n_history.append(model.n.item())  # 记录曼宁系数
        losses.append(loss.item())  # 记录损失

        mu[0] = torch.tensor(S0_history, device=device, dtype=torch.float32).mean()  # 更新S0均值
        mu[1] = torch.tensor(n_history, device=device, dtype=torch.float32).mean()  # 更新n均值
        sigma[0] = max((max(S0_history) - min(S0_history)) / 6, 1e-6)  # 更新S0方差
        sigma[1] = max((max(n_history) - min(n_history)) / 6, 1e-6)  # 更新n方差

        val_mse = compute_validation_mse(model, x_val, t_val, h_val, Q_val)  # 计算验证MSE
        val_mses.append(val_mse)  # 记录验证MSE

        if epoch % 200 == 0:  # 每200次打印结果
            print(f"Epoch {epoch}, Loss: {loss.item():.6f}, Val MSE: {val_mse:.6f}, "
                  f"S0: {model.S0.item():.6f}, n: {model.n.item():.6f}")

    print(f"ADAM Training Time: {time() - start_time:.2f}s")  # 打印ADAM阶段时间

    print("Training with L-BFGS on GPU (5000 epochs)...")
    start_time = time()  # 记录L-BFGS阶段开始时间

    def closure():  # 定义L-BFGS闭包函数
        optimizer_lbfgs.zero_grad()  # 清零梯度
        loss = loss_fn_C(model, x_u, t_u, h_u, Q_u, x_c, t_c, mu, sigma)  # 计算总损失
        if torch.isnan(loss):  # 检查是否出现NaN
            print("NaN detected in L-BFGS")
            return loss
        loss.backward()  # 反向传播
        return loss

    for epoch in range(5000):  # 5000次L-BFGS迭代
        loss = optimizer_lbfgs.step(closure)  # 执行一步优化
        if torch.isnan(loss):  # 检查是否出现NaN
            print(f"NaN detected at epoch {epoch + 2000}")
            break
        S0_history.append(model.S0.item())  # 记录床坡
        n_history.append(model.n.item())  # 记录曼宁系数
        losses.append(loss.item())  # 记录损失

        mu[0] = torch.tensor(S0_history, device=device, dtype=torch.float32).mean()  # 更新S0均值
        mu[1] = torch.tensor(n_history, device=device, dtype=torch.float32).mean()  # 更新n均值
        sigma[0] = max((max(S0_history) - min(S0_history)) / 6, 1e-6)  # 更新S0方差
        sigma[1] = max((max(n_history) - min(n_history)) / 6, 1e-6)  # 更新n方差

        val_mse = compute_validation_mse(model, x_val, t_val, h_val, Q_val)  # 计算验证MSE
        val_mses.append(val_mse)  # 记录验证MSE

        if (epoch + 2000) % 500 == 0:  # 每500次打印结果
            print(f"Epoch {epoch + 2000}, Loss: {loss.item():.6f}, Val MSE: {val_mse:.6f}, "
                  f"S0: {model.S0.item():.6f}, n: {model.n.item():.6f}")

    print(f"L-BFGS Training Time: {time() - start_time:.2f}s")  # 打印L-BFGS阶段时间
    return losses, val_mses, S0_history, n_history  # 返回训练历史


# 可视化结果
def plot_results(model, x_u, t_u, h_u, Q_u, x_val, t_val, losses, val_mses, S0_history, n_history):
    """绘制训练结果"""
    x_u, t_u, h_u = x_u.to('cpu'), t_u.to('cpu'), h_u.to('cpu')  # 数据移到CPU用于绘图
    dh_dt_u = torch.zeros_like(h_u)  # 初始化 Δh_t
    dh_dt_u[1:] = (h_u[1:] - h_u[:-1]) / 10  # 计算数值导数
    dh_dt_u[0] = dh_dt_u[1]  # 边界填充
    with torch.no_grad():  # 无梯度计算
        model = model.to('cpu')  # 模型移到CPU
        h_pred, Q_pred = model(x_u, t_u, dh_dt_u)  # 预测 h 和 Q

    h_u, Q_u = h_u.numpy(), Q_u.numpy()  # 转换为NumPy数组
    h_pred, Q_pred = h_pred.numpy(), Q_pred.numpy()  # 同上
    t_u = t_u.numpy()  # 同上

    plt.figure(figsize=(15, 10))  # 设置图形大小
    plt.subplot(2, 2, 1)  # 子图1:损失曲线
    plt.plot(losses, label='Loss')  # 绘制总损失
    plt.plot(val_mses, label='Validation MSE')  # 绘制验证MSE
    plt.xlabel('Epoch')  # x轴标签
    plt.ylabel('Value')  # y轴标签
    plt.yscale('log')  # 对数刻度
    plt.legend()  # 添加图例

    plt.subplot(2, 2, 2)  # 子图2:流量预测
    plt.plot(t_u[:len(t_u) // 2], Q_u[:len(t_u) // 2], label='True Q (x=0)')  # 真实流量
    plt.plot(t_u[:len(t_u) // 2], Q_pred[:len(t_u) // 2], '--', label='Pred Q (x=0)')  # 预测流量
    plt.xlabel('Time (s)')  # x轴标签
    plt.ylabel('Flow Rate (m^3/s)')  # y轴标签
    plt.legend()  # 添加图例

    plt.subplot(2, 2, 3)  # 子图3:水深预测
    plt.plot(t_u[:len(t_u) // 2], h_u[:len(t_u) // 2], label='True h (x=0)')  # 真实水深
    plt.plot(t_u[:len(t_u) // 2], h_pred[:len(t_u) // 2], '--', label='Pred h (x=0)')  # 预测水深
    plt.xlabel('Time (s)')  # x轴标签
    plt.ylabel('Water Depth (m)')  # y轴标签
    plt.legend()  # 添加图例

    plt.subplot(2, 2, 4)  # 子图4:参数估计
    plt.plot(S0_history, label='S0')  # 绘制床坡估计
    plt.plot(n_history, label='n')  # 绘制曼宁系数估计
    plt.axhline(y=0.0025, color='r', linestyle='--', label='True S0')  # 真实床坡
    plt.axhline(y=0.02, color='g', linestyle='--', label='True n')  # 真实曼宁系数
    plt.xlabel('Epoch')  # x轴标签
    plt.ylabel('Parameter Value')  # y轴标签
    plt.legend()  # 添加图例

    plt.tight_layout()  # 调整布局
    plt.show()  # 显示图形


# 主程序
if __name__ == "__main__":
    (x_u, t_u, h_u, Q_u), (x_c, t_c), (x_val, t_val) = generate_data()  # 生成数据
    h_val = h_u[:300]  # 设置验证集水深
    Q_val = Q_u[:300]  # 设置验证集流量

    model = PINN_C().to(device)  # 初始化模型并移到GPU
    losses, val_mses, S0_history, n_history = train_model(model, x_u, t_u, h_u, Q_u, x_c, t_c, x_val, t_val, h_val,
                                                          Q_val)  # 训练模型

    plot_results(model, x_u, t_u, h_u, Q_u, x_val, t_val, losses, val_mses, S0_history, n_history)  # 可视化结果

    print(f"Final S0: {model.S0.item():.6f}, True S0: 0.0025")  # 打印最终床坡估计
    print(f"Final n: {model.n.item():.6f}, True n: 0.02")  # 打印最终曼宁系数估计
    print(f"Final Validation MSE: {val_mses[-1]:.6f}")  # 打印最终验证MSE

结论

本论文将PINNs作为河流水流率的替代水文模型,研究其在模拟环境中的性能。结果表明PINNs可同化数据并发现控制系统动力学的微分方程参数,能规避概念模型中参数估计的劣势和经验方法中数据缺乏的问题。且在不同条件下,架构A和架构C的PINNs模型分别有较好表现。

后续展望

作者认为未来可将PINNs应用于解决流域的实时系统。此外,进一步探索PINNs在不同真实世界场景中的应用和深入研究物理和数据学习平衡任务的理论。之后可以优化架构C中参数方差的值以提高系统预测能力。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/2304913.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

Python Django系列—入门实例(二)

数据库配置 现在&#xff0c;打开 mysite/settings.py 。这是个包含了 Django 项目设置的 Python 模块。 默认情况下&#xff0c;​ DATABASES 配置使用 SQLite。如果你是数据库新手&#xff0c;或者只是想尝试 Django&#xff0c;这是最简单的选择。SQLite 包含在 Python 中…

【DeepSeek系列】05 DeepSeek核心算法改进点总结

文章目录 一、DeepSeek概要二、4个重要改进点2.1 多头潜在注意力2.2 混合专家模型MoE2.3 多Token预测3.4 GRPO强化学习策略 三、2个重要思考3.1 大规模强化学习3.2 蒸馏方法&#xff1a;小模型也可以很强大 一、DeepSeek概要 2024年&#xff5e;2025年初&#xff0c;DeepSeek …

独立开发者之Google Analytics使用教程

Google Analytics&#xff08;GA&#xff09;是Google提供的一款免费的网络分析服务&#xff0c;用于追踪和报告网站流量。以下是独立开发者如何使用Google Analytics的详细教程&#xff1a; 1. 创建Google Analytics账户 注册Google账户&#xff1a;如果你还没有Google账户&…

C++ 编程语言简介

C 是一种通用编程语言&#xff0c;它是作为 C 语言的增强而开发的&#xff0c;以包含面向对象的范例。它是一种命令式和编译语言。 C 是一种高级的通用编程语言&#xff0c;专为系统和应用程序编程而设计。它由贝尔实验室的 Bjarne Stroustrup 于 1983 年开发&#xff0c;作为…

计算机毕业设计SpringBoot+Vue.js明星周边产品销售网站(源码+文档+PPT+讲解)

温馨提示&#xff1a;文末有 CSDN 平台官方提供的学长联系方式的名片&#xff01; 温馨提示&#xff1a;文末有 CSDN 平台官方提供的学长联系方式的名片&#xff01; 温馨提示&#xff1a;文末有 CSDN 平台官方提供的学长联系方式的名片&#xff01; 作者简介&#xff1a;Java领…

使用Windbg调试目标进程排查C++软件异常的一般步骤与要点分享

目录 1、概述 2、将Windbg附加到已经启动起来的目标进程上&#xff0c;或者用Windbg启动目标程序 2.1、将Windbg附加到已经启动起来的目标进程上 2.2、用Windbg启动目标程序 2.3、Windbg关联到目标进程上会中断下来&#xff0c;输入g命令将该中断跳过去 3、分析实例说明 …

ddd 文章总结分享,ddd实战代码分享, 领域驱动设计java实战源码大全,我看过的ddd java源码

1. 前段时间研究ddd, 收藏了很多相关知识&#xff0c;分享出来&#xff0c;希望能够帮助更多的小伙伴了解ddd, 什么是领域驱动设计&#xff0c;并分享在github发现的开源ddd代码 2. ddd 必须强烈点赞阿里两位大佬&#xff0c;一个为殷浩&#xff0c; 一个为cola作者 2.1.1 殷浩…

什么是MySql的主从复制(主从同步)?

主页还有其他面试题总结&#xff0c;有需要的可以去看一下&#xff0c;喜欢的就留个三连再走吧~ 1.什么是MySql的主从复制原理&#xff1f; 主从复制的核心就是二进制binlog&#xff08;DDL&#xff08;数据定义语言&#xff09;语句和DML&#xff08;数据操纵语言&#xff09…

蓝桥云课python代码

第一章语言基础 第一节编程基础 1 python开发环境 第一个Python程序 # 打印"Hello World" print("Hello World")# 打印2的100次方 print(2 ** 100)# 打印112 print("11",1 1)""" Hello World 126765060022822940149670320537…

c#丰田PLC ToyoPuc TCP协议快速读写 to c# Toyota PLC ToyoPuc读写

源代码下载 <------下载地址 历史背景与发展 TOYOPUC协议源于丰田工机&#xff08;TOYODA&#xff09;的自动化技术积累。丰田工机成立于1941年&#xff0c;最初是丰田汽车的机床部门&#xff0c;后独立为专注于工业机械与控制系统的公司。2006年与光洋精工&#xff08;Ko…

深入解析-无状态服务-StatefulSet (一)

一、有状态服务 VS 无状态服务 1.无状态服务介绍 1.数据方面&#xff1a;无状态服务不会在本地存储持久化数据.多个实例可以共享相同的持久化数据 2.结果方面&#xff1a;多个服务实例对于同一个用户请求的响应结果是完全一致的 3.关系方面&#xff1a;这种多服务实例之间是…

hackmyvm-buster

题目地址 信息收集 主机发现 ┌──(root㉿kali)-[/home/kali] └─# arp-scan -I eth1 192.168.56.0/24 Interface: eth1, type: EN10MB, MAC: 00:0c:29:34:da:f5, IPv4: 192.168.56.103 WARNING: Cannot open MAC/Vendor file ieee-oui.txt: Permission denied WARNING: C…

【原创】Windows11安装WSL“无法解析服务器的名称或地址”问题解决方法

原因分析 出现这个问题一开始以为WSL设置了某个服务器&#xff0c;但是通过运行 nslookup www.microsoft.com 出现下面的提示 PS C:\Windows\system32> nslookup www.microsoft.com 服务器: UnKnown Address: 2408:8000:XXXX:2b00:8:8:8:8非权威应答: 名称: e13678…

网页制作08-html,css,javascript初认识のhtml使用框架结构,请先建立站点!

框架一般由框架集和框架组成。 框架集就像一个大的容器&#xff0c;包括所有的框架&#xff0c;是框架的集合。 框架是框架集中一个独立的区域用于显示一个独立的网页文档。 框架集是文件html&#xff0c;它定义一组框架的布局和属性&#xff0c;包括框架的数目&#xff0c;框架…

【Vscode 使用】集合1

一、使用make工具管理工程 windows下&#xff0c;下载mingw64&#xff0c;配置好mingw64\bin 为 Win10系统全局变量后。 在mingw64/bin目录下找到mingw32-make.exe工具。复制一份改名为&#xff1a;make.exe&#xff0c;没错&#xff0c;就是那么简单&#xff0c;mingw64自带m…

文章精读篇——用于遥感小样本语义分割的可学习Prompt

题目&#xff1a;Learnable Prompt for Few-Shot Semantic Segmentation in Remote Sensing Domain 会议&#xff1a;CVPR 2024 Workshop 论文&#xff1a;10.48550/arXiv.2404.10307 相关竞赛&#xff1a;https://codalab.lisn.upsaclay.fr/competitions/17568 年份&#…

解决 kubeasz 安装k8s集群跨节点pod 无法使用cluster ip通讯问题

问题描述 使用kubeasz搭建k8s集群后使用的配置文件 # etcd cluster should have odd member(s) (1,3,5,...) [etcd] 192.168.xx.22# master node(s) [kube_master] 192.168.xx.22# work node(s) [kube_node] 192.168.xx.9 192.168.xx.22# [optional] harbor server, a privat…

Docker 搭建 Nginx 服务器

系列文章目录 Docker 搭建 Nginx 服务器 系列文章目录前言一、准备工作二、设置 Nginx 容器的目录结构三、启动一个临时的 Nginx 容器来复制配置文件四、复制 Nginx 配置文件到本地目录五、删除临时 Nginx 容器六、创建并运行 Nginx 容器&#xff0c;挂载本地目录七、修改 ngin…

Spring AI + 大模型开发应用

JAVA SpringAI 大模型开发AI应用DEMO 前言JAVA项目创建示例 前言 在当今快速发展的技术领域&#xff0c;人工智能&#xff08;AI&#xff09;已经成为推动创新和变革的重要力量。然而&#xff0c;AI应用的开发过程往往复杂且耗时&#xff0c;需要开发者具备深厚的技术背景和丰…

【C++11】 并发⽀持库

&#x1f308; 个人主页&#xff1a;Zfox_ &#x1f525; 系列专栏&#xff1a;C从入门到精通 目录 前言&#xff1a;&#x1f680; 并发⽀持库一&#xff1a;&#x1f525; thread库 二&#xff1a;&#x1f525; this_thread 三&#xff1a;&#x1f525; mutex 四&#xff1…