跳转至

PhyLSTM

# linux
wget -nc https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyLSTM/data_boucwen.mat
# windows
# curl https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyLSTM/data_boucwen.mat -o data_boucwen.mat
python phylstm2.py
# linux
wget -nc https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyLSTM/data_boucwen.mat
# windows
# curl https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyLSTM/data_boucwen.mat -o data_boucwen.mat
python phylstm3.py
# linux
wget -nc https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyLSTM/data_boucwen.mat
# windows
# curl https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyLSTM/data_boucwen.mat -o data_boucwen.mat
python phylstm2.py mode=eval EVAL.pretrained_model_path=https://paddle-org.bj.bcebos.com/paddlescience/models/phylstm/phylstm2_pretrained.pdparams
# linux
wget -nc https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyLSTM/data_boucwen.mat
# windows
# curl https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyLSTM/data_boucwen.mat -o data_boucwen.mat
python phylstm3.py mode=eval EVAL.pretrained_model_path=https://paddle-org.bj.bcebos.com/paddlescience/models/phylstm/phylstm3_pretrained.pdparams
预训练模型 指标
phylstm2_pretrained.pdparams loss(sup_valid): 0.00799
phylstm3_pretrained.pdparams loss(sup_valid): 0.03098

1. 背景简介

我们引入了一种创新的物理知识LSTM框架,用于对缺乏数据的非线性结构系统进行元建模。基本概念是将可用但尚不完整的物理知识(如物理定律、科学原理)整合到深度长短时记忆(LSTM)网络中,该网络在可行的解决方案空间内限制和促进学习。物理约束嵌入在损失函数中,以强制执行模型训练,即使在可用训练数据集非常有限的情况下,也能准确地捕捉潜在的系统非线性。特别是对于动态结构,考虑运动方程的物理定律、状态依赖性和滞后本构关系来构建物理损失。嵌入式物理可以缓解过拟合问题,减少对大型训练数据集的需求,并提高训练模型的鲁棒性,使其具有外推能力,从而进行更可靠的预测。因此,物理知识指导的深度学习范式优于传统的非物理指导的数据驱动神经网络。

2. 问题定义

结构系统的元建模旨在开发低保真度(或低阶)模型,以有效地捕捉潜在的非线性输入-输出行为。元模型可以在从高保真度模拟或实际系统感知获得的数据集上进行训练。为了更好地说明,我们考虑一个建筑类型结构并假设地震动力学由低保真度非线性运动方程(EOM)支配:

\[ \mathbf{M} \ddot{\mathbf{u}}+\underbrace{\mathbf{C} \dot{\mathbf{u}}+\lambda \mathbf{K u}+(1-\lambda) \mathbf{K r}}_{\mathbf{h}}=-\mathbf{M} \Gamma a_g \]

其中M是质量矩阵;C为阻尼矩阵;K为刚度矩阵。

控制方程可以改写成一个更一般的形式:

\[ \ddot{\mathbf{u}}+\mathrm{g}=-\Gamma a_g \]

3. 问题求解

接下来开始讲解如何将问题一步一步地转化为 PaddleScience 代码,用深度学习的方法求解该问题。 为了快速理解 PaddleScience,接下来仅对模型构建、方程构建、计算域构建等关键步骤进行阐述,而其余细节请参考 API文档

3.1 模型构建

在 PhyLSTM 问题中,建立 LSTM 网络 Deep LSTM network,用 PaddleScience 代码表示如下

model = ppsci.arch.DeepPhyLSTM(
    cfg.MODEL.input_size,
    eta.shape[2],
    cfg.MODEL.hidden_size,
    cfg.MODEL.model_type,
)

DeepPhyLSTM 参数 input_size 是输入大小,output_size 是输出大小,hidden_size 是隐藏层大小,model_type是模型类型。

3.2 数据构建

运行本问题代码前请按照下方命令下载 data_boucwen.mat

wget -nc -P ./ https://paddle-org.bj.bcebos.com/paddlescience/datasets/PhyLSTM/data_boucwen.mat

本案例涉及读取数据构建,如下所示

mat = scipy.io.loadmat(cfg.DATA_FILE_PATH)
ag_data = mat["input_tf"]  # ag, ad, av
u_data = mat["target_X_tf"]
ut_data = mat["target_Xd_tf"]
utt_data = mat["target_Xdd_tf"]
ag_data = ag_data.reshape([ag_data.shape[0], ag_data.shape[1], 1])
u_data = u_data.reshape([u_data.shape[0], u_data.shape[1], 1])
ut_data = ut_data.reshape([ut_data.shape[0], ut_data.shape[1], 1])
utt_data = utt_data.reshape([utt_data.shape[0], utt_data.shape[1], 1])

t = mat["time"]
dt = t[0, 1] - t[0, 0]

ag_all = ag_data
u_all = u_data
u_t_all = ut_data
u_tt_all = utt_data

# finite difference
N = u_data.shape[1]
phi1 = np.concatenate(
    [
        np.array([-3 / 2, 2, -1 / 2]),
        np.zeros([N - 3]),
    ]
)
temp1 = np.concatenate([-1 / 2 * np.identity(N - 2), np.zeros([N - 2, 2])], axis=1)
temp2 = np.concatenate([np.zeros([N - 2, 2]), 1 / 2 * np.identity(N - 2)], axis=1)
phi2 = temp1 + temp2
phi3 = np.concatenate(
    [
        np.zeros([N - 3]),
        np.array([1 / 2, -2, 3 / 2]),
    ]
)
phi_t0 = (
    1
    / dt
    * np.concatenate(
        [
            np.reshape(phi1, [1, phi1.shape[0]]),
            phi2,
            np.reshape(phi3, [1, phi3.shape[0]]),
        ],
        axis=0,
    )
)
phi_t0 = np.reshape(phi_t0, [1, N, N])

ag_star = ag_all[0:10]
eta_star = u_all[0:10]
eta_t_star = u_t_all[0:10]
eta_tt_star = u_tt_all[0:10]
ag_c_star = ag_all[0:50]
lift_star = -ag_c_star

eta = eta_star
ag = ag_star
lift = lift_star
eta_t = eta_t_star
eta_tt = eta_tt_star
ag_c = ag_c_star
g = -eta_tt - ag
phi_t = np.repeat(phi_t0, ag_c_star.shape[0], axis=0)

3.3 约束构建

设置训练数据集和损失计算函数,返回字段,代码如下所示:

sup_constraint_pde = ppsci.constraint.SupervisedConstraint(
    {
        "dataset": {
            "name": "NamedArrayDataset",
            "input": input_dict_train,
            "label": label_dict_train,
        },
        "sampler": {
            "name": "BatchSampler",
            "drop_last": True,
            "shuffle": True,
        },
        "batch_size": 1,
        "num_workers": 0,
    },
    ppsci.loss.FunctionalLoss(functions.train_loss_func2),
    {
        "eta_pred": lambda out: out["eta_pred"],
        "eta_dot_pred": lambda out: out["eta_dot_pred"],
        "g_pred": lambda out: out["g_pred"],
        "eta_t_pred_c": lambda out: out["eta_t_pred_c"],
        "eta_dot_pred_c": lambda out: out["eta_dot_pred_c"],
        "lift_pred_c": lambda out: out["lift_pred_c"],
    },
    name="sup_train",
)
constraint_pde = {sup_constraint_pde.name: sup_constraint_pde}

3.4 评估器构建

设置评估数据集和损失计算函数,返回字段,代码如下所示:

sup_validator_pde = ppsci.validate.SupervisedValidator(
    {
        "dataset": {
            "name": "NamedArrayDataset",
            "input": input_dict_val,
            "label": label_dict_val,
        },
        "sampler": {
            "name": "BatchSampler",
            "drop_last": False,
            "shuffle": False,
        },
        "batch_size": 1,
        "num_workers": 0,
    },
    ppsci.loss.FunctionalLoss(functions.train_loss_func2),
    {
        "eta_pred": lambda out: out["eta_pred"],
        "eta_dot_pred": lambda out: out["eta_dot_pred"],
        "g_pred": lambda out: out["g_pred"],
        "eta_t_pred_c": lambda out: out["eta_t_pred_c"],
        "eta_dot_pred_c": lambda out: out["eta_dot_pred_c"],
        "lift_pred_c": lambda out: out["lift_pred_c"],
    },
    metric={"metric": ppsci.metric.FunctionalMetric(functions.metric_expr)},
    name="sup_valid",
)
validator_pde = {sup_validator_pde.name: sup_validator_pde}

3.5 超参数设定

接下来我们需要指定训练轮数,此处我们按实验经验,使用 100 轮训练轮数。

epochs: 100

3.6 优化器构建

训练过程会调用优化器来更新模型参数,此处选择 Adam 优化器并设定 learning_rate 为 1e-3。

optimizer = ppsci.optimizer.Adam(cfg.TRAIN.learning_rate)(model)

3.7 模型训练与评估

完成上述设置之后,只需要将上述实例化的对象按顺序传递给 ppsci.solver.Solver

solver = ppsci.solver.Solver(
    model,
    constraint_pde,
    cfg.output_dir,
    optimizer,
    None,
    cfg.TRAIN.epochs,
    cfg.TRAIN.iters_per_epoch,
    save_freq=cfg.TRAIN.save_freq,
    log_freq=cfg.log_freq,
    seed=cfg.seed,
    validator=validator_pde,
    checkpoint_path=cfg.TRAIN.checkpoint_path,
    eval_with_no_grad=cfg.EVAL.eval_with_no_grad,
)

最后启动训练、评估即可:

# train model
solver.train()
# evaluate after finished training
solver.eval()

4. 完整代码

phylstm2.py
# Copyright (c) 2023 PaddlePaddle Authors. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
Reference: https://github.com/zhry10/PhyLSTM.git
"""

from os import path as osp

import functions
import hydra
import numpy as np
import scipy.io
from omegaconf import DictConfig

import ppsci
from ppsci.utils import logger


def train(cfg: DictConfig):
    # set random seed for reproducibility
    ppsci.utils.misc.set_random_seed(cfg.seed)
    # initialize logger
    logger.init_logger("ppsci", osp.join(cfg.output_dir, "train.log"), "info")

    mat = scipy.io.loadmat(cfg.DATA_FILE_PATH)
    ag_data = mat["input_tf"]  # ag, ad, av
    u_data = mat["target_X_tf"]
    ut_data = mat["target_Xd_tf"]
    utt_data = mat["target_Xdd_tf"]
    ag_data = ag_data.reshape([ag_data.shape[0], ag_data.shape[1], 1])
    u_data = u_data.reshape([u_data.shape[0], u_data.shape[1], 1])
    ut_data = ut_data.reshape([ut_data.shape[0], ut_data.shape[1], 1])
    utt_data = utt_data.reshape([utt_data.shape[0], utt_data.shape[1], 1])

    t = mat["time"]
    dt = t[0, 1] - t[0, 0]

    ag_all = ag_data
    u_all = u_data
    u_t_all = ut_data
    u_tt_all = utt_data

    # finite difference
    N = u_data.shape[1]
    phi1 = np.concatenate(
        [
            np.array([-3 / 2, 2, -1 / 2]),
            np.zeros([N - 3]),
        ]
    )
    temp1 = np.concatenate([-1 / 2 * np.identity(N - 2), np.zeros([N - 2, 2])], axis=1)
    temp2 = np.concatenate([np.zeros([N - 2, 2]), 1 / 2 * np.identity(N - 2)], axis=1)
    phi2 = temp1 + temp2
    phi3 = np.concatenate(
        [
            np.zeros([N - 3]),
            np.array([1 / 2, -2, 3 / 2]),
        ]
    )
    phi_t0 = (
        1
        / dt
        * np.concatenate(
            [
                np.reshape(phi1, [1, phi1.shape[0]]),
                phi2,
                np.reshape(phi3, [1, phi3.shape[0]]),
            ],
            axis=0,
        )
    )
    phi_t0 = np.reshape(phi_t0, [1, N, N])

    ag_star = ag_all[0:10]
    eta_star = u_all[0:10]
    eta_t_star = u_t_all[0:10]
    eta_tt_star = u_tt_all[0:10]
    ag_c_star = ag_all[0:50]
    lift_star = -ag_c_star

    eta = eta_star
    ag = ag_star
    lift = lift_star
    eta_t = eta_t_star
    eta_tt = eta_tt_star
    ag_c = ag_c_star
    g = -eta_tt - ag
    phi_t = np.repeat(phi_t0, ag_c_star.shape[0], axis=0)

    model = ppsci.arch.DeepPhyLSTM(
        cfg.MODEL.input_size,
        eta.shape[2],
        cfg.MODEL.hidden_size,
        cfg.MODEL.model_type,
    )
    model.register_input_transform(functions.transform_in)
    model.register_output_transform(functions.transform_out)

    dataset_obj = functions.Dataset(eta, eta_t, g, ag, ag_c, lift, phi_t)
    (
        input_dict_train,
        label_dict_train,
        input_dict_val,
        label_dict_val,
    ) = dataset_obj.get(cfg.TRAIN.epochs)

    sup_constraint_pde = ppsci.constraint.SupervisedConstraint(
        {
            "dataset": {
                "name": "NamedArrayDataset",
                "input": input_dict_train,
                "label": label_dict_train,
            },
            "sampler": {
                "name": "BatchSampler",
                "drop_last": True,
                "shuffle": True,
            },
            "batch_size": 1,
            "num_workers": 0,
        },
        ppsci.loss.FunctionalLoss(functions.train_loss_func2),
        {
            "eta_pred": lambda out: out["eta_pred"],
            "eta_dot_pred": lambda out: out["eta_dot_pred"],
            "g_pred": lambda out: out["g_pred"],
            "eta_t_pred_c": lambda out: out["eta_t_pred_c"],
            "eta_dot_pred_c": lambda out: out["eta_dot_pred_c"],
            "lift_pred_c": lambda out: out["lift_pred_c"],
        },
        name="sup_train",
    )
    constraint_pde = {sup_constraint_pde.name: sup_constraint_pde}

    sup_validator_pde = ppsci.validate.SupervisedValidator(
        {
            "dataset": {
                "name": "NamedArrayDataset",
                "input": input_dict_val,
                "label": label_dict_val,
            },
            "sampler": {
                "name": "BatchSampler",
                "drop_last": False,
                "shuffle": False,
            },
            "batch_size": 1,
            "num_workers": 0,
        },
        ppsci.loss.FunctionalLoss(functions.train_loss_func2),
        {
            "eta_pred": lambda out: out["eta_pred"],
            "eta_dot_pred": lambda out: out["eta_dot_pred"],
            "g_pred": lambda out: out["g_pred"],
            "eta_t_pred_c": lambda out: out["eta_t_pred_c"],
            "eta_dot_pred_c": lambda out: out["eta_dot_pred_c"],
            "lift_pred_c": lambda out: out["lift_pred_c"],
        },
        metric={"metric": ppsci.metric.FunctionalMetric(functions.metric_expr)},
        name="sup_valid",
    )
    validator_pde = {sup_validator_pde.name: sup_validator_pde}

    # initialize solver
    optimizer = ppsci.optimizer.Adam(cfg.TRAIN.learning_rate)(model)
    solver = ppsci.solver.Solver(
        model,
        constraint_pde,
        cfg.output_dir,
        optimizer,
        None,
        cfg.TRAIN.epochs,
        cfg.TRAIN.iters_per_epoch,
        save_freq=cfg.TRAIN.save_freq,
        log_freq=cfg.log_freq,
        seed=cfg.seed,
        validator=validator_pde,
        checkpoint_path=cfg.TRAIN.checkpoint_path,
        eval_with_no_grad=cfg.EVAL.eval_with_no_grad,
    )

    # train model
    solver.train()
    # evaluate after finished training
    solver.eval()


def evaluate(cfg: DictConfig):
    # set random seed for reproducibility
    ppsci.utils.misc.set_random_seed(cfg.seed)
    # initialize logger
    logger.init_logger("ppsci", osp.join(cfg.output_dir, "eval.log"), "info")

    mat = scipy.io.loadmat(cfg.DATA_FILE_PATH)
    ag_data = mat["input_tf"]  # ag, ad, av
    u_data = mat["target_X_tf"]
    ut_data = mat["target_Xd_tf"]
    utt_data = mat["target_Xdd_tf"]
    ag_data = ag_data.reshape([ag_data.shape[0], ag_data.shape[1], 1])
    u_data = u_data.reshape([u_data.shape[0], u_data.shape[1], 1])
    ut_data = ut_data.reshape([ut_data.shape[0], ut_data.shape[1], 1])
    utt_data = utt_data.reshape([utt_data.shape[0], utt_data.shape[1], 1])

    t = mat["time"]
    dt = t[0, 1] - t[0, 0]

    ag_all = ag_data
    u_all = u_data
    u_t_all = ut_data
    u_tt_all = utt_data

    # finite difference
    N = u_data.shape[1]
    phi1 = np.concatenate(
        [
            np.array([-3 / 2, 2, -1 / 2]),
            np.zeros([N - 3]),
        ]
    )
    temp1 = np.concatenate([-1 / 2 * np.identity(N - 2), np.zeros([N - 2, 2])], axis=1)
    temp2 = np.concatenate([np.zeros([N - 2, 2]), 1 / 2 * np.identity(N - 2)], axis=1)
    phi2 = temp1 + temp2
    phi3 = np.concatenate(
        [
            np.zeros([N - 3]),
            np.array([1 / 2, -2, 3 / 2]),
        ]
    )
    phi_t0 = (
        1
        / dt
        * np.concatenate(
            [
                np.reshape(phi1, [1, phi1.shape[0]]),
                phi2,
                np.reshape(phi3, [1, phi3.shape[0]]),
            ],
            axis=0,
        )
    )
    phi_t0 = np.reshape(phi_t0, [1, N, N])

    ag_star = ag_all[0:10]
    eta_star = u_all[0:10]
    eta_t_star = u_t_all[0:10]
    eta_tt_star = u_tt_all[0:10]
    ag_c_star = ag_all[0:50]
    lift_star = -ag_c_star

    eta = eta_star
    ag = ag_star
    lift = lift_star
    eta_t = eta_t_star
    eta_tt = eta_tt_star
    ag_c = ag_c_star
    g = -eta_tt - ag
    phi_t = np.repeat(phi_t0, ag_c_star.shape[0], axis=0)

    model = ppsci.arch.DeepPhyLSTM(
        cfg.MODEL.input_size,
        eta.shape[2],
        cfg.MODEL.hidden_size,
        cfg.MODEL.model_type,
    )
    model.register_input_transform(functions.transform_in)
    model.register_output_transform(functions.transform_out)

    dataset_obj = functions.Dataset(eta, eta_t, g, ag, ag_c, lift, phi_t)
    (
        _,
        _,
        input_dict_val,
        label_dict_val,
    ) = dataset_obj.get(1)

    sup_validator_pde = ppsci.validate.SupervisedValidator(
        {
            "dataset": {
                "name": "NamedArrayDataset",
                "input": input_dict_val,
                "label": label_dict_val,
            },
            "sampler": {
                "name": "BatchSampler",
                "drop_last": False,
                "shuffle": False,
            },
            "batch_size": 1,
            "num_workers": 0,
        },
        ppsci.loss.FunctionalLoss(functions.train_loss_func2),
        {
            "eta_pred": lambda out: out["eta_pred"],
            "eta_dot_pred": lambda out: out["eta_dot_pred"],
            "g_pred": lambda out: out["g_pred"],
            "eta_t_pred_c": lambda out: out["eta_t_pred_c"],
            "eta_dot_pred_c": lambda out: out["eta_dot_pred_c"],
            "lift_pred_c": lambda out: out["lift_pred_c"],
        },
        metric={"metric": ppsci.metric.FunctionalMetric(functions.metric_expr)},
        name="sup_valid",
    )
    validator_pde = {sup_validator_pde.name: sup_validator_pde}

    # initialize solver
    solver = ppsci.solver.Solver(
        model,
        output_dir=cfg.output_dir,
        seed=cfg.seed,
        validator=validator_pde,
        pretrained_model_path=cfg.EVAL.pretrained_model_path,
        eval_with_no_grad=cfg.EVAL.eval_with_no_grad,
    )
    # evaluate
    solver.eval()


@hydra.main(version_base=None, config_path="./conf", config_name="phylstm2.yaml")
def main(cfg: DictConfig):
    if cfg.mode == "train":
        train(cfg)
    elif cfg.mode == "eval":
        evaluate(cfg)
    else:
        raise ValueError(f"cfg.mode should in ['train', 'eval'], but got '{cfg.mode}'")


if __name__ == "__main__":
    main()
phylstm3.py
# Copyright (c) 2023 PaddlePaddle Authors. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
Reference: https://github.com/zhry10/PhyLSTM.git
"""

from os import path as osp

import functions
import hydra
import numpy as np
import scipy.io
from omegaconf import DictConfig

import ppsci
from ppsci.utils import logger


def train(cfg: DictConfig):
    # set random seed for reproducibility
    ppsci.utils.misc.set_random_seed(cfg.seed)
    # initialize logger
    logger.init_logger("ppsci", osp.join(cfg.output_dir, "train.log"), "info")

    mat = scipy.io.loadmat(cfg.DATA_FILE_PATH)
    t = mat["time"]
    dt = 0.02
    n1 = int(dt / 0.005)
    t = t[::n1]

    ag_data = mat["input_tf"][:, ::n1]  # ag, ad, av
    u_data = mat["target_X_tf"][:, ::n1]
    ut_data = mat["target_Xd_tf"][:, ::n1]
    utt_data = mat["target_Xdd_tf"][:, ::n1]
    ag_data = ag_data.reshape([ag_data.shape[0], ag_data.shape[1], 1])
    u_data = u_data.reshape([u_data.shape[0], u_data.shape[1], 1])
    ut_data = ut_data.reshape([ut_data.shape[0], ut_data.shape[1], 1])
    utt_data = utt_data.reshape([utt_data.shape[0], utt_data.shape[1], 1])

    ag_pred = mat["input_pred_tf"][:, ::n1]  # ag, ad, av
    u_pred = mat["target_pred_X_tf"][:, ::n1]
    ut_pred = mat["target_pred_Xd_tf"][:, ::n1]
    utt_pred = mat["target_pred_Xdd_tf"][:, ::n1]
    ag_pred = ag_pred.reshape([ag_pred.shape[0], ag_pred.shape[1], 1])
    u_pred = u_pred.reshape([u_pred.shape[0], u_pred.shape[1], 1])
    ut_pred = ut_pred.reshape([ut_pred.shape[0], ut_pred.shape[1], 1])
    utt_pred = utt_pred.reshape([utt_pred.shape[0], utt_pred.shape[1], 1])

    N = u_data.shape[1]
    phi1 = np.concatenate(
        [
            np.array([-3 / 2, 2, -1 / 2]),
            np.zeros([N - 3]),
        ]
    )
    temp1 = np.concatenate([-1 / 2 * np.identity(N - 2), np.zeros([N - 2, 2])], axis=1)
    temp2 = np.concatenate([np.zeros([N - 2, 2]), 1 / 2 * np.identity(N - 2)], axis=1)
    phi2 = temp1 + temp2
    phi3 = np.concatenate(
        [
            np.zeros([N - 3]),
            np.array([1 / 2, -2, 3 / 2]),
        ]
    )
    phi_t0 = (
        1
        / dt
        * np.concatenate(
            [
                np.reshape(phi1, [1, phi1.shape[0]]),
                phi2,
                np.reshape(phi3, [1, phi3.shape[0]]),
            ],
            axis=0,
        )
    )
    phi_t0 = np.reshape(phi_t0, [1, N, N])

    ag_star = ag_data
    eta_star = u_data
    eta_t_star = ut_data
    eta_tt_star = utt_data
    ag_c_star = np.concatenate([ag_data, ag_pred[0:53]])
    lift_star = -ag_c_star

    eta = eta_star
    ag = ag_star
    lift = lift_star
    eta_t = eta_t_star
    eta_tt = eta_tt_star
    g = -eta_tt - ag
    ag_c = ag_c_star

    phi_t = np.repeat(phi_t0, ag_c_star.shape[0], axis=0)

    model = ppsci.arch.DeepPhyLSTM(
        cfg.MODEL.input_size,
        eta.shape[2],
        cfg.MODEL.hidden_size,
        cfg.MODEL.model_type,
    )
    model.register_input_transform(functions.transform_in)
    model.register_output_transform(functions.transform_out)

    dataset_obj = functions.Dataset(eta, eta_t, g, ag, ag_c, lift, phi_t)
    (
        input_dict_train,
        label_dict_train,
        input_dict_val,
        label_dict_val,
    ) = dataset_obj.get(cfg.TRAIN.epochs)

    sup_constraint_pde = ppsci.constraint.SupervisedConstraint(
        {
            "dataset": {
                "name": "NamedArrayDataset",
                "input": input_dict_train,
                "label": label_dict_train,
            },
            "sampler": {
                "name": "BatchSampler",
                "drop_last": False,
                "shuffle": False,
            },
            "batch_size": 1,
            "num_workers": 0,
        },
        ppsci.loss.FunctionalLoss(functions.train_loss_func3),
        {
            "eta_pred": lambda out: out["eta_pred"],
            "eta_dot_pred": lambda out: out["eta_dot_pred"],
            "g_pred": lambda out: out["g_pred"],
            "eta_t_pred_c": lambda out: out["eta_t_pred_c"],
            "eta_dot_pred_c": lambda out: out["eta_dot_pred_c"],
            "lift_pred_c": lambda out: out["lift_pred_c"],
            "g_t_pred_c": lambda out: out["g_t_pred_c"],
            "g_dot_pred_c": lambda out: out["g_dot_pred_c"],
        },
        name="sup_train",
    )
    constraint_pde = {sup_constraint_pde.name: sup_constraint_pde}

    sup_validator_pde = ppsci.validate.SupervisedValidator(
        {
            "dataset": {
                "name": "NamedArrayDataset",
                "input": input_dict_val,
                "label": label_dict_val,
            },
            "sampler": {
                "name": "BatchSampler",
                "drop_last": False,
                "shuffle": False,
            },
            "batch_size": 1,
            "num_workers": 0,
        },
        ppsci.loss.FunctionalLoss(functions.train_loss_func3),
        {
            "eta_pred": lambda out: out["eta_pred"],
            "eta_dot_pred": lambda out: out["eta_dot_pred"],
            "g_pred": lambda out: out["g_pred"],
            "eta_t_pred_c": lambda out: out["eta_t_pred_c"],
            "eta_dot_pred_c": lambda out: out["eta_dot_pred_c"],
            "lift_pred_c": lambda out: out["lift_pred_c"],
            "g_t_pred_c": lambda out: out["g_t_pred_c"],
            "g_dot_pred_c": lambda out: out["g_dot_pred_c"],
        },
        metric={"metric": ppsci.metric.FunctionalMetric(functions.metric_expr)},
        name="sup_valid",
    )
    validator_pde = {sup_validator_pde.name: sup_validator_pde}

    # initialize solver
    optimizer = ppsci.optimizer.Adam(cfg.TRAIN.learning_rate)(model)
    solver = ppsci.solver.Solver(
        model,
        constraint_pde,
        cfg.output_dir,
        optimizer,
        None,
        cfg.TRAIN.epochs,
        cfg.TRAIN.iters_per_epoch,
        save_freq=cfg.TRAIN.save_freq,
        log_freq=cfg.log_freq,
        seed=cfg.seed,
        validator=validator_pde,
        checkpoint_path=cfg.TRAIN.checkpoint_path,
        eval_with_no_grad=cfg.EVAL.eval_with_no_grad,
    )

    # train model
    solver.train()
    # evaluate after finished training
    solver.eval()


def evaluate(cfg: DictConfig):
    # set random seed for reproducibility
    ppsci.utils.misc.set_random_seed(cfg.seed)
    # initialize logger
    logger.init_logger("ppsci", osp.join(cfg.output_dir, "train.log"), "info")

    mat = scipy.io.loadmat(cfg.DATA_FILE_PATH)
    t = mat["time"]
    dt = 0.02
    n1 = int(dt / 0.005)
    t = t[::n1]

    ag_data = mat["input_tf"][:, ::n1]  # ag, ad, av
    u_data = mat["target_X_tf"][:, ::n1]
    ut_data = mat["target_Xd_tf"][:, ::n1]
    utt_data = mat["target_Xdd_tf"][:, ::n1]
    ag_data = ag_data.reshape([ag_data.shape[0], ag_data.shape[1], 1])
    u_data = u_data.reshape([u_data.shape[0], u_data.shape[1], 1])
    ut_data = ut_data.reshape([ut_data.shape[0], ut_data.shape[1], 1])
    utt_data = utt_data.reshape([utt_data.shape[0], utt_data.shape[1], 1])

    ag_pred = mat["input_pred_tf"][:, ::n1]  # ag, ad, av
    u_pred = mat["target_pred_X_tf"][:, ::n1]
    ut_pred = mat["target_pred_Xd_tf"][:, ::n1]
    utt_pred = mat["target_pred_Xdd_tf"][:, ::n1]
    ag_pred = ag_pred.reshape([ag_pred.shape[0], ag_pred.shape[1], 1])
    u_pred = u_pred.reshape([u_pred.shape[0], u_pred.shape[1], 1])
    ut_pred = ut_pred.reshape([ut_pred.shape[0], ut_pred.shape[1], 1])
    utt_pred = utt_pred.reshape([utt_pred.shape[0], utt_pred.shape[1], 1])

    N = u_data.shape[1]
    phi1 = np.concatenate(
        [
            np.array([-3 / 2, 2, -1 / 2]),
            np.zeros([N - 3]),
        ]
    )
    temp1 = np.concatenate([-1 / 2 * np.identity(N - 2), np.zeros([N - 2, 2])], axis=1)
    temp2 = np.concatenate([np.zeros([N - 2, 2]), 1 / 2 * np.identity(N - 2)], axis=1)
    phi2 = temp1 + temp2
    phi3 = np.concatenate(
        [
            np.zeros([N - 3]),
            np.array([1 / 2, -2, 3 / 2]),
        ]
    )
    phi_t0 = (
        1
        / dt
        * np.concatenate(
            [
                np.reshape(phi1, [1, phi1.shape[0]]),
                phi2,
                np.reshape(phi3, [1, phi3.shape[0]]),
            ],
            axis=0,
        )
    )
    phi_t0 = np.reshape(phi_t0, [1, N, N])

    ag_star = ag_data
    eta_star = u_data
    eta_t_star = ut_data
    eta_tt_star = utt_data
    ag_c_star = np.concatenate([ag_data, ag_pred[0:53]])
    lift_star = -ag_c_star

    eta = eta_star
    ag = ag_star
    lift = lift_star
    eta_t = eta_t_star
    eta_tt = eta_tt_star
    g = -eta_tt - ag
    ag_c = ag_c_star

    phi_t = np.repeat(phi_t0, ag_c_star.shape[0], axis=0)

    model = ppsci.arch.DeepPhyLSTM(
        cfg.MODEL.input_size,
        eta.shape[2],
        cfg.MODEL.hidden_size,
        cfg.MODEL.model_type,
    )
    model.register_input_transform(functions.transform_in)
    model.register_output_transform(functions.transform_out)

    dataset_obj = functions.Dataset(eta, eta_t, g, ag, ag_c, lift, phi_t)
    (
        _,
        _,
        input_dict_val,
        label_dict_val,
    ) = dataset_obj.get(1)

    sup_validator_pde = ppsci.validate.SupervisedValidator(
        {
            "dataset": {
                "name": "NamedArrayDataset",
                "input": input_dict_val,
                "label": label_dict_val,
            },
            "sampler": {
                "name": "BatchSampler",
                "drop_last": False,
                "shuffle": False,
            },
            "batch_size": 1,
            "num_workers": 0,
        },
        ppsci.loss.FunctionalLoss(functions.train_loss_func3),
        {
            "eta_pred": lambda out: out["eta_pred"],
            "eta_dot_pred": lambda out: out["eta_dot_pred"],
            "g_pred": lambda out: out["g_pred"],
            "eta_t_pred_c": lambda out: out["eta_t_pred_c"],
            "eta_dot_pred_c": lambda out: out["eta_dot_pred_c"],
            "lift_pred_c": lambda out: out["lift_pred_c"],
            "g_t_pred_c": lambda out: out["g_t_pred_c"],
            "g_dot_pred_c": lambda out: out["g_dot_pred_c"],
        },
        metric={"metric": ppsci.metric.FunctionalMetric(functions.metric_expr)},
        name="sup_valid",
    )
    validator_pde = {sup_validator_pde.name: sup_validator_pde}

    # initialize solver
    solver = ppsci.solver.Solver(
        model,
        output_dir=cfg.output_dir,
        seed=cfg.seed,
        validator=validator_pde,
        pretrained_model_path=cfg.EVAL.pretrained_model_path,
        eval_with_no_grad=cfg.EVAL.eval_with_no_grad,
    )

    # evaluate
    solver.eval()


@hydra.main(version_base=None, config_path="./conf", config_name="phylstm3.yaml")
def main(cfg: DictConfig):
    if cfg.mode == "train":
        train(cfg)
    elif cfg.mode == "eval":
        evaluate(cfg)
    else:
        raise ValueError(f"cfg.mode should in ['train', 'eval'], but got '{cfg.mode}'")


if __name__ == "__main__":
    main()

5. 结果展示

PhyLSTM2 案例针对 epoch=100 和 learning_rate=1e-3 的参数配置进行了实验,结果返回Loss为 0.00799。

PhyLSTM3 案例针对 epoch=200 和 learning_rate=1e-3 的参数配置进行了实验,结果返回Loss为 0.03098。

6. 参考资料