主题043：机器学习与数据驱动仿真

1. 引言

1.1 传统仿真的挑战

传统的多物理场耦合仿真方法虽然精确，但面临以下挑战：

计算成本高昂

• 复杂几何的高分辨率网格划分
• 长时间瞬态分析需要大量时间步
• 参数研究和优化设计需要多次仿真
• 实时仿真难以实现

模型复杂度限制

• 微观尺度现象难以直接建模
• 多尺度耦合计算量巨大
• 不确定性量化需要大量样本
• 高维参数空间探索困难

实验数据利用不足

• 传统方法难以直接融合实验数据
• 模型修正依赖人工经验
• 数据驱动发现物理规律的能力有限
  

1.2 机器学习的机遇

机器学习为多物理场仿真带来新的可能性：

加速计算

• 代理模型替代昂贵仿真
• 降阶模型减少自由度
• 实时预测和快速优化
• 不确定性传播的高效计算

数据融合

• 实验与仿真数据联合学习
• 自动模型修正和校准
• 缺失数据的智能填补
• 多源异构数据整合

新发现能力

• 从数据中发现控制方程
• 识别隐藏物理规律
• 构建数据驱动的本构模型
• 发现新的材料行为

1.3 本章内容概述

本章系统介绍机器学习在多物理场仿真中的应用，包括：

1. 机器学习基础：监督学习、无监督学习、强化学习
2. 神经网络理论：前馈网络、卷积网络、循环网络
3. 物理信息神经网络：将物理约束嵌入神经网络
4. 降阶模型：本征正交分解、动态模态分解
5. 数据驱动仿真：从数据学习物理系统
6. 工程应用：热传导、流体力学、结构力学等

---

2. 机器学习基础理论

2.1 机器学习分类

2.1.1 监督学习

监督学习从标注数据中学习输入到输出的映射关系。

回归问题

给定训练数据集 $\{({\mathbf{x}}_i,y_i){\}}_{i=1}^N$，学习函数 $f:\mathcal{X}\to \mathcal{Y}$：

$$y=f(\mathbf{x};\mathbf{\theta })+ϵ$$

其中 $\mathbf{\theta }$ 为模型参数，$ϵ$ 为噪声。

常用损失函数：

• 均方误差（MSE）：$\mathcal{L}=\frac{1}{N}{\sum }_{i=1}^N(y_i-{\hat{y}}_i{)}^2$
• 平均绝对误差（MAE）：$\mathcal{L}=\frac{1}{N}{\sum }_{i=1}^N|y_i-{\hat{y}}_i|$
• Huber损失：结合MSE和MAE的优点

分类问题

输出为离散类别标签，常用交叉熵损失：

$$\mathcal{L}=-\frac{1}{N}\sum _{i=1}^N\sum _{c=1}^C{y}_{i,c}\log ({\hat{y}}_{i,c})$$

2.1.2 无监督学习

无监督学习从无标注数据中发现隐藏结构。

聚类分析

将数据划分为若干组，使得组内相似度高，组间相似度低：

$$J=\sum _{k=1}^K\sum _{\mathbf{x}\in C_k}\Vert \mathbf{x}-{\mathbf{\mu }}_k{\Vert }^2$$

降维技术

主成分分析（PCA）寻找数据的主要变化方向：

$$\mathbf{z}={\mathbf{W}}^T\mathbf{x}$$

其中 $\mathbf{W}$ 包含协方差矩阵的前 $k$ 个特征向量。

2.1.3 强化学习

强化学习通过与环境交互学习最优策略。

马尔可夫决策过程

状态转移：$s_{t+1}\sim P(\cdot |s_t,a_t)$

奖励函数：$r_t=R(s_t,a_t)$

目标：最大化累积奖励 $\mathbb{E}[{\sum }_{t=0}^T{\gamma }^t{r}_t]$

2.2 经典机器学习算法

2.2.1 线性回归与多项式回归

线性回归

模型：$y={\mathbf{w}}^T\mathbf{x}+b$

解析解：$\mathbf{w}=({\mathbf{X}}^T\mathbf{X}{)}^{-1}{\mathbf{X}}^T\mathbf{y}$

多项式回归

通过特征映射实现非线性拟合：

$$\varphi (\mathbf{x})=[1,x_1,x_2,x_1^2,x_1{x}_2,x_2^2,...]$$

2.2.2 支持向量机

分类SVM

寻找最优超平面：${\mathbf{w}}^T\mathbf{x}+b=0$

优化问题：

$$\underset{\mathbf{w},b}{\min }\frac{1}{2}\Vert \mathbf{w}{\Vert }^2+C\sum _{i=1}^N\max (0,1-y_i({\mathbf{w}}^T{\mathbf{x}}_i+b))$$

核技巧

通过核函数隐式映射到高维空间：

$$K({\mathbf{x}}_i,{\mathbf{x}}_j)=\varphi ({\mathbf{x}}_i{)}^T\varphi ({\mathbf{x}}_j)$$

常用核函数：

• 线性核：$K(\mathbf{x},\mathbf{y})={\mathbf{x}}^T\mathbf{y}$
• 多项式核：$K(\mathbf{x},\mathbf{y})=(\gamma {\mathbf{x}}^T\mathbf{y}+r{)}^d$
• RBF核：$K(\mathbf{x},\mathbf{y})=\exp (-\gamma \Vert \mathbf{x}-\mathbf{y}{\Vert }^2)$

2.2.3 决策树与集成方法

决策树

通过递归划分特征空间构建树形结构。

划分准则：

• 信息增益（分类）：$IG=H(D)-{\sum }_{v\in Values(A)}\frac{|D_v|}{|D|}H(D_v)$
• 基尼指数：$Gini(D)=1-{\sum }_{k=1}^K{p}_k^2$
• 方差减少（回归）

随机森林

集成多棵决策树，通过Bagging降低方差：

$$\hat{y}=\frac{1}{M}\sum _{m=1}^M{T}_m(\mathbf{x})$$

梯度提升树

通过Boosting顺序训练弱学习器：

$$F_m(\mathbf{x})=F_{m-1}(\mathbf{x})+\nu \cdot h_m(\mathbf{x})$$

2.2.4 高斯过程回归

高斯过程定义

函数值的有限集合服从联合高斯分布：

$$\mathbf{f}\sim \mathcal{N}(0,\mathbf{K})$$

其中 $K_{ij}=k({\mathbf{x}}_i,{\mathbf{x}}_j)$ 为核函数。

预测分布

给定训练数据，新点的预测：

$$p(f_{\ast }|\mathbf{X},\mathbf{y},{\mathbf{x}}_{\ast })=\mathcal{N}({\mu }_{\ast },{\sigma }_{\ast }^2)$$

其中：

$${\mu }_{\ast }={\mathbf{k}}_{\ast }^T{\mathbf{K}}^{-1}\mathbf{y}$$
$${\sigma }_{\ast }^2=k({\mathbf{x}}_{\ast },{\mathbf{x}}_{\ast })-{\mathbf{k}}_{\ast }^T{\mathbf{K}}^{-1}{\mathbf{k}}_{\ast }$$

优势：提供预测不确定性量化

2.3 模型评估与选择

2.3.1 交叉验证

K折交叉验证

1. 将数据分为K个子集
2. 每次使用K-1个子集训练，1个子集验证
3. 重复K次，确保每个子集都作为验证集
4. 平均K次验证结果

留一法交叉验证

K等于样本数，适用于小数据集。

2.3.2 性能指标

回归指标

• 均方根误差（RMSE）：$RMSE=\sqrt{\frac{1}{N}{\sum }_{i=1}^N(y_i-{\hat{y}}_i{)}^2}$
• 决定系数（$R^2$）：$R^2=1-\frac{\sum (y_i-{\hat{y}}_i{)}^2}{\sum (y_i-\bar{y}{)}^2}$
• 平均绝对百分比误差（MAPE）

分类指标

• 准确率：$Accuracy=\frac{TP+TN}{TP+TN+FP+FN}$
• 精确率：$Precision=\frac{TP}{TP+FP}$
• 召回率：$Recall=\frac{TP}{TP+FN}$
• F1分数：$F1=2\cdot \frac{Precision\cdot Recall}{Precision+Recall}$

2.3.3 超参数优化

网格搜索

在预定义的参数网格上穷举搜索：

$${\mathbf{\theta }}^{\ast }=\arg \underset{\mathbf{\theta }\in \Theta }{\min }CV(\mathbf{\theta })$$

随机搜索

在参数空间中随机采样，更高效：

$${\mathbf{\theta }}_i\sim p(\mathbf{\theta }),i=1,...,n$$

贝叶斯优化

构建目标函数的概率代理模型，选择最有希望的点：

$${\mathbf{x}}_{next}=\arg \underset{\mathbf{x}}{\max }\alpha (\mathbf{x};\mathcal{D})$$

其中 $\alpha$ 为采集函数（Expected Improvement, Upper Confidence Bound等）。

---

3. 神经网络与深度学习

3.1 前馈神经网络

3.1.1 网络结构

基本单元

神经元计算：$z={\sum }_{j=1}^d{w}_j{x}_j+b={\mathbf{w}}^T\mathbf{x}+b$

激活函数：$a=\sigma (z)$

多层感知机（MLP）

$${\mathbf{h}}^{(l)}=\sigma ({\mathbf{W}}^{(l)}{\mathbf{h}}^{(l-1)}+{\mathbf{b}}^{(l)})$$

$$\hat{y}={\mathbf{W}}^{(L)}{\mathbf{h}}^{(L-1)}+{\mathbf{b}}^{(L)}$$

通用近似定理

具有足够多隐藏层神经元的前馈神经网络可以以任意精度近似任意连续函数。

3.1.2 激活函数

Sigmoid

$$\sigma (x)=\frac{1}{1+e^{-x}}$$

特点：输出范围(0,1)，适合二分类；但存在梯度消失问题。

Tanh

$$\tanh (x)=\frac{e^x-e^{-x}}{e^x+e^{-x}}$$

特点：输出范围(-1,1)，零中心化；仍有梯度消失。

ReLU（Rectified Linear Unit）

$$\text{ReLU}(x)=\max (0,x)$$

特点：计算简单，缓解梯度消失；但存在"死亡ReLU"问题。

Leaky ReLU

$$\text{LeakyReLU}(x)=\max (\alpha x,x),\alpha >0$$

特点：解决死亡ReLU问题。

Swish

$$\text{Swish}(x)=x\cdot \sigma (x)=\frac{x}{1+e^{-x}}$$

特点：自门控激活函数，性能优异。

3.1.3 损失函数与优化

常用损失函数

• 回归：MSE, MAE, Huber Loss
• 分类：交叉熵, Focal Loss
• 物理约束：PDE残差, 边界条件损失

梯度下降

批量梯度下降：

$${\mathbf{\theta }}_{t+1}={\mathbf{\theta }}_t-\eta {\nabla }_{\mathbf{\theta }}\mathcal{L}({\mathbf{\theta }}_t)$$

随机梯度下降（SGD）：

$${\mathbf{\theta }}_{t+1}={\mathbf{\theta }}_t-\eta {\nabla }_{\mathbf{\theta }}\mathcal{L}({\mathbf{\theta }}_t;{\mathbf{x}}_i,y_i)$$

小批量梯度下降：

$${\mathbf{\theta }}_{t+1}={\mathbf{\theta }}_t-\eta \frac{1}{m}\sum _{i=1}^m{\nabla }_{\mathbf{\theta }}\mathcal{L}({\mathbf{\theta }}_t;{\mathbf{x}}_i,y_i)$$

自适应优化算法

Adam（Adaptive Moment Estimation）

一阶矩估计：${\mathbf{m}}_t={\beta }_1{\mathbf{m}}_{t-1}+(1-{\beta }_1){\mathbf{g}}_t$

二阶矩估计：${\mathbf{v}}_t={\beta }_2{\mathbf{v}}_{t-1}+(1-{\beta }_2){\mathbf{g}}_t^2$

偏差修正：

$${\hat{\mathbf{m}}}_t=\frac{{\mathbf{m}}_t}{1-{\beta }_1^t},{\hat{\mathbf{v}}}_t=\frac{{\mathbf{v}}_t}{1-{\beta }_2^t}$$

参数更新：

$${\mathbf{\theta }}_{t+1}={\mathbf{\theta }}_t-\frac{\eta }{\sqrt{{\hat{\mathbf{v}}}_t}+ϵ}{\hat{\mathbf{m}}}_t$$

3.1.4 正则化技术

L1/L2正则化

L2正则化（权重衰减）：${\mathcal{L}}_{reg}=\mathcal{L}+\lambda \Vert \mathbf{W}{\Vert }_F^2$

L1正则化（稀疏性）：${\mathcal{L}}_{reg}=\mathcal{L}+\lambda \Vert \mathbf{W}{\Vert }_1$

Dropout

训练时随机丢弃神经元：

$${\mathbf{h}}^{(l)}={\mathbf{m}}^{(l)}\odot \sigma ({\mathbf{W}}^{(l)}{\mathbf{h}}^{(l-1)}+{\mathbf{b}}^{(l)})$$

其中 $m_j^{(l)}\sim \text{Bernoulli}(p)$

早停（Early Stopping）

监控验证集性能，在过拟合开始前停止训练。

批归一化（Batch Normalization）

对每层输入进行归一化：

$${\hat{x}}_i=\frac{x_i-{\mu }_B}{\sqrt{{\sigma }_B^2+ϵ}}$$

$$y_i=\gamma {\hat{x}}_i+\beta$$

3.2 卷积神经网络

3.2.1 卷积操作

一维卷积

$$(f\ast g)[n]=\sum _{m=-M}^{M}f[m]g[n-m]$$

二维卷积

$$(I\ast K)[i,j]=\sum _{m=0}^{M-1}\sum _{n=0}^{N-1}I[i+m,j+n]K[m,n]$$

卷积层特性

• 局部连接：每个神经元只连接局部区域
• 权重共享：同一卷积核在整个输入上滑动
• 平移等变性：目标位置变化，特征响应随之移动

3.2.2 池化操作

最大池化

$$y_{i,j}=\underset{(m,n)\in R_{i,j}}{\max }{x}_{m,n}$$

平均池化

$$y_{i,j}=\frac{1}{|R_{i,j}|}\sum _{(m,n)\in R_{i,j}}{x}_{m,n}$$

作用：降低特征图尺寸，增强平移不变性，减少计算量。

3.2.3 经典网络架构

LeNet（1998）

结构：Conv → Pool → Conv → Pool → FC → FC → Output

AlexNet（2012）

创新：ReLU激活，Dropout正则化，GPU并行训练

VGGNet（2014）

特点：使用小卷积核（3×3）堆叠，网络深度达16-19层

ResNet（2015）

残差连接：y=F(x,{Wi})+x\mathbf{y} = \mathcal{F}(\mathbf{x}, \{W_i\}) + \mathbf{x}y=F(x,{Wi​})+x

解决深层网络的梯度消失问题，可训练152+层网络。

3.3 循环神经网络

3.3.1 基本RNN

前向传播

$${\mathbf{h}}_t=\sigma ({\mathbf{W}}_{hh}{\mathbf{h}}_{t-1}+{\mathbf{W}}_{xh}{\mathbf{x}}_t+{\mathbf{b}}_h)$$

$${\mathbf{y}}_t={\mathbf{W}}_{hy}{\mathbf{h}}_t+{\mathbf{b}}_y$$

问题：长期依赖问题，梯度消失/爆炸

3.3.2 LSTM（长短期记忆网络）

门控机制

遗忘门：${\mathbf{f}}_t=\sigma ({\mathbf{W}}_f[{\mathbf{h}}_{t-1},{\mathbf{x}}_t]+{\mathbf{b}}_f)$

输入门：${\mathbf{i}}_t=\sigma ({\mathbf{W}}_i[{\mathbf{h}}_{t-1},{\mathbf{x}}_t]+{\mathbf{b}}_i)$

候选状态：${\tilde{\mathbf{C}}}_t=\tanh ({\mathbf{W}}_C[{\mathbf{h}}_{t-1},{\mathbf{x}}_t]+{\mathbf{b}}_C)$

细胞状态更新：${\mathbf{C}}_t={\mathbf{f}}_t\odot {\mathbf{C}}_{t-1}+{\mathbf{i}}_t\odot {\tilde{\mathbf{C}}}_t$

输出门：${\mathbf{o}}_t=\sigma ({\mathbf{W}}_o[{\mathbf{h}}_{t-1},{\mathbf{x}}_t]+{\mathbf{b}}_o)$

隐藏状态：${\mathbf{h}}_t={\mathbf{o}}_t\odot \tanh ({\mathbf{C}}_t)$

3.3.3 GRU（门控循环单元）

简化版LSTM，合并细胞状态和隐藏状态：

更新门：${\mathbf{z}}_t=\sigma ({\mathbf{W}}_z[{\mathbf{h}}_{t-1},{\mathbf{x}}_t])$

重置门：${\mathbf{r}}_t=\sigma ({\mathbf{W}}_r[{\mathbf{h}}_{t-1},{\mathbf{x}}_t])$

候选状态：${\tilde{\mathbf{h}}}_t=\tanh (\mathbf{W}[{\mathbf{r}}_t\odot {\mathbf{h}}_{t-1},{\mathbf{x}}_t])$

隐藏状态：${\mathbf{h}}_t=(1-{\mathbf{z}}_t)\odot {\mathbf{h}}_{t-1}+{\mathbf{z}}_t\odot {\tilde{\mathbf{h}}}_t$

3.4 自动编码器与生成模型

3.4.1 自动编码器

结构

编码器：$\mathbf{z}=f_{enc}(\mathbf{x})$

解码器：$\hat{\mathbf{x}}=f_{dec}(\mathbf{z})$

损失函数：$\mathcal{L}=\Vert \mathbf{x}-\hat{\mathbf{x}}{\Vert }^2$

应用：降维、特征学习、去噪、异常检测

3.4.2 变分自动编码器（VAE）

概率框架

编码器学习后验分布：$q_{\varphi }(\mathbf{z}|\mathbf{x})$

解码器学习似然：$p_{\theta }(\mathbf{x}|\mathbf{z})$

目标函数（ELBO）：

$$\mathcal{L}(\theta ,\varphi ;\mathbf{x})={\mathbb{E}}_{q_{\varphi }(\mathbf{z}|\mathbf{x})}[\log p_{\theta }(\mathbf{x}|\mathbf{z})]-D_{KL}(q_{\varphi }(\mathbf{z}|\mathbf{x})\Vert p(\mathbf{z}))$$

3.4.3 生成对抗网络（GAN）

对抗训练

生成器：$G(\mathbf{z})$ 从噪声生成样本

判别器：$D(\mathbf{x})$ 区分真实和生成样本

目标函数：

$$\underset{G}{\min }\underset{D}{\max }V(D,G)={\mathbb{E}}_{\mathbf{x}\sim p_{data}}[\log D(\mathbf{x})]+{\mathbb{E}}_{\mathbf{z}\sim p_z}[\log (1-D(G(\mathbf{z})))]$$

---

4. 物理信息神经网络(PINN)

4.1 PINN基本原理

4.1.1 核心思想

物理信息神经网络将物理定律（偏微分方程）作为约束嵌入神经网络训练过程中。

网络表示

神经网络 $u_{NN}(\mathbf{x},t;\mathbf{\theta })$ 近似PDE解 $u(\mathbf{x},t)$

损失函数组成

$$\mathcal{L}={\mathcal{L}}_{data}+{\mathcal{L}}_{PDE}+{\mathcal{L}}_{BC}+{\mathcal{L}}_{IC}$$

其中：

• ${\mathcal{L}}_{data}$：数据拟合损失
• ${\mathcal{L}}_{PDE}$：PDE残差损失
• ${\mathcal{L}}_{BC}$：边界条件损失
• ${\mathcal{L}}_{IC}$：初始条件损失

4.1.2 PDE残差计算

考虑一般形式的PDE：

$$\mathcal{F}[u]=f(\mathbf{x},t),\mathbf{x}\in \Omega ,t\in [0,T]$$

PDE残差：

$$r(\mathbf{x},t;\mathbf{\theta })=\mathcal{F}[u_{NN}(\mathbf{x},t;\mathbf{\theta })]-f(\mathbf{x},t)$$

残差损失：

$${\mathcal{L}}_{PDE}=\frac{1}{N_f}\sum _{i=1}^{N_f}|r({\mathbf{x}}_f^i,t_f^i;\mathbf{\theta }){|}^2$$

自动微分

使用自动微分计算高阶导数：

$$\frac{\partial u_{NN}}{\partial x}=\text{autodiff}(u_{NN},x)$$

$$\frac{{\partial }^2{u}_{NN}}{\partial x^2}=\text{autodiff}(\frac{\partial u_{NN}}{\partial x},x)$$

4.2 PINN求解各类PDE

4.2.1 泊松方程

问题描述

$${\nabla }^{2}u=f(\mathbf{x}),\mathbf{x}\in \Omega$$

边界条件：$u=g(\mathbf{x}),\mathbf{x}\in \partial \Omega$

PINN实现

残差：$r={\nabla }^2{u}_{NN}-f$

边界损失：${\mathcal{L}}_{BC}=\frac{1}{N_b}{\sum }_{i=1}^{N_b}|u_{NN}({\mathbf{x}}_b^i)-g({\mathbf{x}}_b^i){|}^2$

4.2.2 热传导方程

问题描述

$$\frac{\partial u}{\partial t}=\alpha {\nabla }^{2}u+f(\mathbf{x},t)$$

初始条件：$u(\mathbf{x},0)=u_0(\mathbf{x})$

边界条件：$u=g(\mathbf{x},t)$ 或 $\frac{\partial u}{\partial n}=q(\mathbf{x},t)$

PINN实现

残差：$r=\frac{\partial u_{NN}}{\partial t}-\alpha {\nabla }^2{u}_{NN}-f$

初始损失：${\mathcal{L}}_{IC}=\frac{1}{N_{ic}}{\sum }_{i=1}^{N_{ic}}|u_{NN}({\mathbf{x}}_{ic}^i,0)-u_0({\mathbf{x}}_{ic}^i){|}^2$

4.2.3 波动方程

问题描述

$$\frac{{\partial }^{2}u}{\partial t^2}=c^2{\nabla }^{2}u$$

初始条件：$u(\mathbf{x},0)=u_0(\mathbf{x})$, $\frac{\partial u}{\partial t}(\mathbf{x},0)=v_0(\mathbf{x})$

PINN实现

残差：$r=\frac{{\partial }^2{u}_{NN}}{\partial t^2}-c^2{\nabla }^2{u}_{NN}$

4.2.4 纳维-斯托克斯方程

问题描述

连续性方程：$\nabla \cdot \mathbf{u}=0$

动量方程：$\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \nabla )\mathbf{u}=-\frac{1}{\rho }\nabla p+\nu {\nabla }^2\mathbf{u}$

PINN实现

残差1：$r_1=\nabla \cdot {\mathbf{u}}_{NN}$

残差2：$r_2=\frac{\partial {\mathbf{u}}_{NN}}{\partial t}+({\mathbf{u}}_{NN}\cdot \nabla ){\mathbf{u}}_{NN}+\frac{1}{\rho }\nabla p_{NN}-\nu {\nabla }^2{\mathbf{u}}_{NN}$

4.3 PINN的改进与扩展

4.3.1 自适应权重

损失平衡问题

不同损失项的量级差异导致训练困难。

自适应权重策略

1. 梯度统计法：根据梯度大小调整权重
2. 学习率退火：动态调整各损失项权重
3. 神经切线核（NTK）分析：根据收敛速度调整权重

4.3.2 域分解PINN

思想

将计算域分解为多个子域，每个子域使用独立的神经网络。

界面条件

连续性条件：$u_i=u_j$ 在界面 ${\Gamma }_{ij}$ 上

通量连续性：$\frac{\partial u_i}{\partial n}=\frac{\partial u_j}{\partial n}$ 在界面 ${\Gamma }_{ij}$ 上

优势：

• 处理复杂几何
• 并行计算
• 局部细化

4.3.3 时空自适应采样

残差自适应细化（RAR）

1. 在粗网格上训练初始PINN
2. 计算残差分布
3. 在残差大的区域增加配点
4. 重新训练

基于误差的自适应

$${\mathbf{x}}_{new}=\arg \underset{\mathbf{x}}{\max }|r(\mathbf{x};\mathbf{\theta })|$$

4.3.4 参数推断PINN

逆问题求解

同时学习解和未知参数：

$$\underset{\mathbf{\theta },\mathbf{\lambda }}{\min }{\mathcal{L}}_{data}(\mathbf{\theta })+{\mathcal{L}}_{PDE}(\mathbf{\theta },\mathbf{\lambda })$$

其中 $\mathbf{\lambda }$ 为待推断的物理参数（材料属性、源项等）。

应用：

• 材料参数识别
• 源项反演
• 边界条件识别

4.4 PINN的优势与局限

4.4.1 优势

1. 网格无关：不需要传统数值方法的网格划分
2. 高维问题：可处理高维PDE（维数灾难缓解）
3. 逆问题自然处理：可同时求解正问题和逆问题
4. 数据融合：可结合物理定律和实验数据
5. 连续解：解在整个域上连续可微

4.4.2 局限

1. 训练成本：需要大量计算资源训练神经网络
2. 收敛困难：复杂PDE可能难以收敛
3. 边界层问题：陡峭梯度区域精度不足
4. 长时间积分：瞬态问题长时间模拟困难
5. 超参数敏感：网络架构和训练参数需要仔细调优

---

5. 降阶模型与代理模型

5.1 本征正交分解（POD）

5.1.1 POD基本原理

目标

寻找最优基函数，使得降阶表示的误差最小。

快照矩阵

收集 $N_s$ 个快照（解样本）：$\mathbf{S}=[{\mathbf{u}}_1,{\mathbf{u}}_2,...,{\mathbf{u}}_{N_s}]\in {\mathbb{R}}^{N\times N_s}$

奇异值分解

$$\mathbf{S}=\mathbf{U\Sigma }{\mathbf{V}}^T$$

其中 $\mathbf{U}=[{\mathbf{\varphi }}_1,{\mathbf{\varphi }}_2,...,{\mathbf{\varphi }}_N]$ 为POD模态。

降阶近似

$$\mathbf{u}\approx \sum _{i=1}^r{a}_i{\mathbf{\varphi }}_i={\mathbf{\Phi }}_r\mathbf{a}$$

其中 $r\ll N$ 为降阶维度。

能量捕获

前 $r$ 个模态捕获的能量：

$$E_r=\frac{{\sum }_{i=1}^r{\sigma }_i^2}{{\sum }_{i=1}^{\min (N,N_s)}{\sigma }_i^2}$$

5.1.2 Galerkin投影

降阶动力学

将全阶模型投影到POD子空间：

$${\mathbf{\Phi }}_r^T\mathbf{M}{\mathbf{\Phi }}_r\dot{\mathbf{a}}={\mathbf{\Phi }}_r^T\mathbf{f}({\mathbf{\Phi }}_r\mathbf{a})$$

$$\tilde{\mathbf{M}}\dot{\mathbf{a}}=\tilde{\mathbf{f}}(\mathbf{a})$$

离线-在线分解

• 离线阶段：计算POD基，投影算子
• 在线阶段：求解降阶系统（快速）

5.2 动态模态分解（DMD）

5.2.1 DMD算法

基本假设

动力学可由线性算子近似：${\mathbf{u}}_{k+1}=\mathbf{A}{\mathbf{u}}_k$

算法步骤

1. 构建数据矩阵：$\mathbf{X}=[{\mathbf{u}}_1,...,{\mathbf{u}}_{m-1}]$, ${\mathbf{X}}^{\prime }=[{\mathbf{u}}_2,...,{\mathbf{u}}_m]$

2. 近似线性算子：${\mathbf{X}}^{\prime }\approx \mathbf{AX}$

3. 低秩近似：$\tilde{\mathbf{A}}={\mathbf{U}}_r^T{\mathbf{X}}^{\prime }{\mathbf{V}}_r{\mathbf{\Sigma }}_r^{-1}$

4. 特征分解：$\tilde{\mathbf{A}}{\mathbf{w}}_j={\lambda }_j{\mathbf{w}}_j$

5. DMD模态：${\mathbf{\varphi }}_j={\mathbf{X}}^{\prime }{\mathbf{V}}_r{\mathbf{\Sigma }}_r^{-1}{\mathbf{w}}_j$

重构与预测

$$\mathbf{u}(t)=\sum _{j=1}^r{b}_j{\mathbf{\varphi }}_j{e}^{{\omega }_{j}t}$$

其中 ${\omega }_j=\ln ({\lambda }_j)/\Delta t$

5.2.2 扩展DMD方法

精确DMD

直接使用 ${\mathbf{X}}^{\prime }$ 而非 $\mathbf{X}$ 计算模态。

压缩DMD

使用随机投影降低计算成本。

前向后向DMD

结合前向和后向时间演化提高稳定性。

5.3 高斯过程回归代理模型

5.3.1 作为代理模型

应用场景

• 昂贵仿真的快速近似
• 优化设计空间探索
• 不确定性量化

训练与预测

1. 选择训练点：$\mathbf{X}=[{\mathbf{x}}_1,...,{\mathbf{x}}_N]$
2. 运行高保真仿真获得响应：$\mathbf{y}=[y_1,...,y_N]$
3. 训练GP模型
4. 在新点预测：$\mu ({\mathbf{x}}_{\ast })$ 和 $\sigma ({\mathbf{x}}_{\ast })$

自适应采样

基于预测不确定性选择新训练点：

$${\mathbf{x}}_{new}=\arg \underset{\mathbf{x}}{\max }\sigma (\mathbf{x})$$

5.4 神经网络代理模型

5.4.1 全连接网络代理

输入输出映射

输入：设计参数 $\mathbf{x}\in {\mathbb{R}}^d$

输出：系统响应 $\mathbf{y}\in {\mathbb{R}}^m$

网络：$\mathbf{y}=f_{NN}(\mathbf{x};\mathbf{\theta })$

训练数据

通过实验设计（DOE）生成训练样本：

• 拉丁超立方采样（LHS）
• Sobol序列
• 正交数组

5.4.2 卷积神经网络代理

适用场景

输入输出为空间场（图像、网格）。

编码器-解码器架构

编码器：将输入场压缩为低维表示

瓶颈层：参数化表示

解码器：重构输出场

U-Net架构

跳跃连接保留细节信息，适合场到场的映射。

5.5 多保真度建模

5.5.1 多保真度框架

保真度层次

• 低保真（LF）：快速但精度低（简化模型、粗网格）
• 中保真（MF）：平衡精度和速度
• 高保真（HF）：精确但昂贵（精细模型）

协同Kriging

$$y_{HF}(\mathbf{x})=\rho y_{LF}(\mathbf{x})+\delta (\mathbf{x})$$

其中 $\rho$ 为比例系数，$\delta$ 为差异项。

多保真度神经网络

使用低保真数据预训练，高保真数据微调。

---

6. 数据驱动的多物理场仿真

6.1 数据驱动发现控制方程

6.1.1 SINDy算法

稀疏识别非线性动力学（Sparse Identification of Nonlinear Dynamics）

给定数据 {x(ti)}i=1m\{\mathbf{x}(t_i)\}_{i=1}^m{x(ti​)}i=1m​，识别：

$$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x})$$

算法步骤

1. 计算导数：$\dot{\mathbf{X}}=[\dot{\mathbf{x}}(t_1),...,\dot{\mathbf{x}}(t_m){]}^T$

2. 构建库矩阵：$\mathbf{\Theta }(\mathbf{X})=[1,\mathbf{X},{\mathbf{X}}^2,...,\sin (\mathbf{X}),...]$

3. 稀疏回归：$\dot{\mathbf{X}}=\mathbf{\Theta }(\mathbf{X})\mathbf{\Xi }$

4. 使用LASSO求解稀疏系数：${\min }_{\mathbf{\Xi }}\Vert \dot{\mathbf{X}}-\mathbf{\Theta }(\mathbf{X})\mathbf{\Xi }{\Vert }_2^2+\lambda \Vert \mathbf{\Xi }{\Vert }_1$

6.1.2 PDE-FIND

PDE功能识别（PDE Functional Identification of Nonlinear Dynamics）

扩展SINDy到PDE发现：

$$u_t=\mathbf{\Theta }(u,u_x,u_{xx},...)\mathbf{\xi }$$

识别空间导数项的稀疏组合。

6.2 数据驱动的本构建模

6.2.1 神经网络本构模型

传统本构模型

$$\mathbf{\sigma }=\mathbf{f}(\mathbf{\epsilon },{\mathbf{\epsilon }}^p,\alpha ,...)$$

神经网络替代

使用神经网络学习本构关系：

$$\mathbf{\sigma }=NN(\mathbf{\epsilon },{\mathbf{\epsilon }}^p,\alpha ,...;\mathbf{\theta })$$

约束嵌入

• 热力学一致性约束
• 客观性（框架无关性）
• 材料对称性

6.2.2 循环神经网络本构模型

历史依赖性

材料响应依赖于加载历史：

$${\mathbf{\sigma }}_t=RNN({\mathbf{\epsilon }}_t,{\mathbf{h}}_{t-1};\mathbf{\theta })$$

应用

• 粘弹性材料
• 塑性变形
• 损伤演化

6.3 数字孪生中的机器学习

6.3.1 实时状态估计

卡尔曼滤波与神经网络结合

使用神经网络学习卡尔曼增益：

$${\mathbf{K}}_t=NN({\mathbf{P}}_t,{\mathbf{H}}_t;\mathbf{\theta })$$

粒子滤波与深度学习

使用自编码器降维，在隐空间进行粒子滤波。

6.3.2 预测性维护

剩余使用寿命预测

$$RU{L}_t=LSTM({\mathbf{x}}_{t-T:t};\mathbf{\theta })$$

输入：传感器历史数据
输出：剩余寿命预测

异常检测

使用自编码器重构误差检测异常：

$$e=\Vert \mathbf{x}-\hat{\mathbf{x}}\Vert >{ϵ}_{threshold}$$

6.4 多物理场数据融合

6.4.1 多模态学习

数据类型

• 数值仿真数据
• 实验测量数据
• 图像/视频数据
• 文本描述

融合策略

• 早期融合：特征级联
• 中期融合：共享表示
• 晚期融合：决策集成

6.4.2 图神经网络多物理场

图表示

将多物理场系统表示为图：

• 节点：物理量（温度、位移、压力）
• 边：物理耦合关系

消息传递

$${\mathbf{h}}_i^{(l+1)}=\sigma (\sum _{j\in \mathcal{N}(i)}{\mathbf{W}}^{(l)}{\mathbf{h}}_j^{(l)}+{\mathbf{b}}^{(l)})$$

---

7. 工程应用案例

7.1 热传导问题

7.1.1 PINN求解瞬态热传导

问题描述

二维瞬态热传导：$\frac{\partial T}{\partial t}=\alpha {\nabla }^{2}T+Q(\mathbf{x},t)$

边界条件：混合边界（Dirichlet + Neumann）

PINN设置

• 网络：5层全连接，每层50神经元
• 激活函数：tanh
• 训练点：内部10000点，边界2000点
• 优化器：Adam + L-BFGS

结果

与有限元解对比，相对误差 < 1%

7.1.2 热传导代理模型

应用场景

快速预测不同材料参数和边界条件下的温度场。

方法

• 使用POD降维
• 高斯过程回归映射参数到模态系数
• 在线预测速度提升1000倍

7.2 流体力学问题

7.2.1 圆柱绕流PINN

问题描述

不可压缩Navier-Stokes方程，雷诺数Re=100

挑战

• 涡脱落周期性
• 压力-速度耦合
• 长时间稳定性

解决方案

• 分阶段训练
• 周期性边界条件处理
• 压力泊松方程约束

7.2.2 湍流模型

数据驱动湍流模型

使用DNS数据训练RANS湍流模型：

$${\tau }_{ij}^{RANS}=NN(k,\epsilon ,\frac{\partial U_i}{\partial x_j};\mathbf{\theta })$$

7.3 结构力学问题

7.3.1 非线性材料识别

问题

从载荷-位移曲线识别材料参数。

方法

• PINN同时学习位移场和材料参数
• 损失函数包含测量数据拟合和平衡方程

7.3.2 结构损伤检测

方法

• 使用卷积神经网络分析振动信号
• 自编码器学习健康状态特征
• 重构误差指示损伤位置和程度

7.4 多物理场耦合问题

7.4.1 热-力耦合PINN

控制方程

热传导：$\rho c_p\dot{T}=k{\nabla }^{2}T+\beta T_0{\dot{\epsilon }}_{kk}^e$

力学平衡：$\nabla \cdot \mathbf{\sigma }=0$

本构：$\mathbf{\sigma }=\mathbf{C}:(\mathbf{\epsilon }-{\mathbf{\epsilon }}^{th})$

耦合PINN架构

• 两个神经网络：$T_{NN}$ 和 ${\mathbf{u}}_{NN}$
• 共享部分层学习耦合关系
• 联合损失函数包含两个物理场的约束

7.4.2 流-固耦合代理模型

应用场景

心脏瓣膜血流动力学快速预测。

方法

• 使用3D CNN处理时空数据
• 编码器提取流场特征
• 解码器预测瓣膜变形

---

8. Python实现与案例分析

8.1 案例1：PINN求解一维热传导方程

import numpy as np
import torch
import torch.nn as nn
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

class HeatPINN(nn.Module):
    """物理信息神经网络求解1D热传导方程"""
    def __init__(self, layers):
        super(HeatPINN, self).__init__()
        self.layers = nn.ModuleList()
        for i in range(len(layers)-1):
            self.layers.append(nn.Linear(layers[i], layers[i+1]))
        self.activation = nn.Tanh()
    
    def forward(self, x, t):
        """前向传播"""
        X = torch.cat([x, t], dim=1)
        for i, layer in enumerate(self.layers[:-1]):
            X = self.activation(layer(X))
        return self.layers[-1](X)
    
    def net_u(self, x, t):
        """网络输出（温度）"""
        return self.forward(x, t)
    
    def net_f(self, x, t, alpha):
        """计算PDE残差"""
        u = self.net_u(x, t)
        
        # 自动微分计算导数
        u_t = torch.autograd.grad(u, t, grad_outputs=torch.ones_like(u),
                                   create_graph=True)[0]
        u_x = torch.autograd.grad(u, x, grad_outputs=torch.ones_like(u),
                                   create_graph=True)[0]
        u_xx = torch.autograd.grad(u_x, x, grad_outputs=torch.ones_like(u_x),
                                    create_graph=True)[0]
        
        # PDE残差: u_t - alpha * u_xx = 0
        f = u_t - alpha * u_xx
        return f

# 训练函数
def train_pinn(model, x_f, t_f, x_ic, t_ic, u_ic, x_bc, t_bc, u_bc, 
               alpha, epochs=10000, lr=0.001):
    """训练PINN"""
    optimizer = torch.optim.Adam(model.parameters(), lr=lr)
    
    losses = []
    for epoch in range(epochs):
        optimizer.zero_grad()
        
        # PDE残差损失
        f_pred = model.net_f(x_f, t_f, alpha)
        loss_f = torch.mean(f_pred**2)
        
        # 初始条件损失
        u_ic_pred = model.net_u(x_ic, t_ic)
        loss_ic = torch.mean((u_ic_pred - u_ic)**2)
        
        # 边界条件损失
        u_bc_pred = model.net_u(x_bc, t_bc)
        loss_bc = torch.mean((u_bc_pred - u_bc)**2)
        
        # 总损失
        loss = loss_f + loss_ic + loss_bc
        
        loss.backward()
        optimizer.step()
        
        losses.append(loss.item())
        
        if epoch % 1000 == 0:
            print(f'Epoch {epoch}, Loss: {loss.item():.6f}, '
                  f'PDE: {loss_f.item():.6f}, IC: {loss_ic.item():.6f}, BC: {loss_bc.item():.6f}')
    
    return losses

# 主程序
def case1_heat_pinn():
    """案例1：PINN求解1D热传导方程"""
    print("="*60)
    print("案例1：PINN求解一维热传导方程")
    print("="*60)
    
    # 问题参数
    L = 1.0  # 杆长
    T = 0.5  # 总时间
    alpha = 0.1  # 热扩散系数
    
    # 生成训练数据
    np.random.seed(42)
    
    # 内部配点（PDE残差计算点）
    N_f = 10000
    x_f = np.random.rand(N_f, 1) * L
    t_f = np.random.rand(N_f, 1) * T
    
    # 初始条件点
    N_ic = 200
    x_ic = np.linspace(0, L, N_ic).reshape(-1, 1)
    t_ic = np.zeros_like(x_ic)
    u_ic = np.sin(np.pi * x_ic / L)  # 初始温度分布
    
    # 边界条件点
    N_bc = 200
    t_bc = np.linspace(0, T, N_bc).reshape(-1, 1)
    x_bc_left = np.zeros_like(t_bc)
    x_bc_right = np.ones_like(t_bc) * L
    u_bc_left = np.zeros_like(t_bc)
    u_bc_right = np.zeros_like(t_bc)
    
    # 转换为torch张量
    x_f = torch.tensor(x_f, dtype=torch.float32, requires_grad=True)
    t_f = torch.tensor(t_f, dtype=torch.float32, requires_grad=True)
    x_ic = torch.tensor(x_ic, dtype=torch.float32)
    t_ic = torch.tensor(t_ic, dtype=torch.float32)
    u_ic = torch.tensor(u_ic, dtype=torch.float32)
    x_bc = torch.tensor(np.vstack([x_bc_left, x_bc_right]), dtype=torch.float32)
    t_bc = torch.tensor(np.vstack([t_bc, t_bc]), dtype=torch.float32)
    u_bc = torch.tensor(np.vstack([u_bc_left, u_bc_right]), dtype=torch.float32)
    
    # 创建PINN模型
    layers = [2, 50, 50, 50, 50, 1]
    model = HeatPINN(layers)
    
    print("\n开始训练PINN...")
    losses = train_pinn(model, x_f, t_f, x_ic, t_ic, u_ic, x_bc, t_bc, u_bc,
                        alpha, epochs=10000, lr=0.001)
    
    # 生成预测网格
    nx, nt = 100, 100
    x_grid = np.linspace(0, L, nx)
    t_grid = np.linspace(0, T, nt)
    X, T_grid = np.meshgrid(x_grid, t_grid)
    
    x_test = torch.tensor(X.flatten(), dtype=torch.float32).reshape(-1, 1)
    t_test = torch.tensor(T_grid.flatten(), dtype=torch.float32).reshape(-1, 1)
    
    # 预测
    model.eval()
    with torch.no_grad():
        u_pred = model.net_u(x_test, t_test).numpy().reshape(nt, nx)
    
    # 解析解（对比）
    u_exact = np.sin(np.pi * X / L) * np.exp(-alpha * (np.pi/L)**2 * T_grid)
    
    # 计算误差
    error = np.abs(u_pred - u_exact)
    print(f"\n最大误差: {np.max(error):.6f}")
    print(f"平均误差: {np.mean(error):.6f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    # 训练损失
    axes[0, 0].semilogy(losses)
    axes[0, 0].set_xlabel('Epoch')
    axes[0, 0].set_ylabel('Loss')
    axes[0, 0].set_title('Training Loss')
    axes[0, 0].grid(True)
    
    # PINN预测
    im1 = axes[0, 1].contourf(X, T_grid, u_pred, levels=50, cmap='hot')
    axes[0, 1].set_xlabel('x')
    axes[0, 1].set_ylabel('t')
    axes[0, 1].set_title('PINN Prediction')
    plt.colorbar(im1, ax=axes[0, 1])
    
    # 解析解
    im2 = axes[0, 2].contourf(X, T_grid, u_exact, levels=50, cmap='hot')
    axes[0, 2].set_xlabel('x')
    axes[0, 2].set_ylabel('t')
    axes[0, 2].set_title('Exact Solution')
    plt.colorbar(im2, ax=axes[0, 2])
    
    # 误差分布
    im3 = axes[1, 0].contourf(X, T_grid, error, levels=50, cmap='coolwarm')
    axes[1, 0].set_xlabel('x')
    axes[1, 0].set_ylabel('t')
    axes[1, 0].set_title('Absolute Error')
    plt.colorbar(im3, ax=axes[1, 0])
    
    # 特定时刻的温度分布
    time_indices = [0, nt//4, nt//2, 3*nt//4, nt-1]
    for idx in time_indices:
        axes[1, 1].plot(x_grid, u_pred[idx, :], label=f't={t_grid[idx]:.2f}')
    axes[1, 1].set_xlabel('x')
    axes[1, 1].set_ylabel('Temperature')
    axes[1, 1].set_title('Temperature Profiles at Different Times')
    axes[1, 1].legend()
    axes[1, 1].grid(True)
    
    # 中心点温度随时间变化
    center_idx = nx // 2
    axes[1, 2].plot(t_grid, u_pred[:, center_idx], 'b-', label='PINN')
    axes[1, 2].plot(t_grid, u_exact[:, center_idx], 'r--', label='Exact')
    axes[1, 2].set_xlabel('t')
    axes[1, 2].set_ylabel('Temperature')
    axes[1, 2].set_title('Temperature at Center Point')
    axes[1, 2].legend()
    axes[1, 2].grid(True)
    
    plt.tight_layout()
    plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题043\\output\\case1_heat_pinn.png', dpi=150)
    plt.close()
    
    print("\n案例1完成！")
    return model, u_pred, u_exact

# 运行案例1
if __name__ == "__main__":
    case1_heat_pinn()

8.2 案例2：POD降阶模型

def case2_pod_rom():
    """案例2：本征正交分解降阶模型"""
    print("\n" + "="*60)
    print("案例2：POD降阶模型")
    print("="*60)
    
    # 生成高保真数据（2D热传导）
    nx, ny = 50, 50
    nt = 100
    Lx, Ly = 1.0, 1.0
    T_final = 1.0
    alpha = 0.01
    
    dx, dy = Lx/(nx-1), Ly/(ny-1)
    dt = T_final / nt
    
    # 空间网格
    x = np.linspace(0, Lx, nx)
    y = np.linspace(0, Ly, ny)
    X, Y = np.meshgrid(x, y)
    
    # 参数变化：不同边界条件
    n_samples = 50
    snapshots = []
    params = []
    
    print("\n生成高保真快照数据...")
    for i in range(n_samples):
        # 随机边界条件参数
        T_left = np.random.uniform(0.5, 1.5)
        T_right = np.random.uniform(0.5, 1.5)
        T_top = np.random.uniform(0.5, 1.5)
        T_bottom = np.random.uniform(0.5, 1.5)
        params.append([T_left, T_right, T_top, T_bottom])
        
        # 稳态热传导求解（有限差分）
        T = np.zeros((ny, nx))
        T[:, 0] = T_left
        T[:, -1] = T_right
        T[0, :] = T_bottom
        T[-1, :] = T_top
        
        # 迭代求解
        for _ in range(1000):
            T_old = T.copy()
            T[1:-1, 1:-1] = 0.25 * (T_old[1:-1, :-2] + T_old[1:-1, 2:] + 
                                     T_old[:-2, 1:-1] + T_old[2:, 1:-1])
            # 保持边界条件
            T[:, 0] = T_left
            T[:, -1] = T_right
            T[0, :] = T_bottom
            T[-1, :] = T_top
        
        snapshots.append(T.flatten())
    
    snapshots = np.array(snapshots)  # (n_samples, nx*ny)
    params = np.array(params)
    
    print(f"快照矩阵形状: {snapshots.shape}")
    
    # POD分析
    print("\n执行POD分析...")
    
    # 计算均值并去均值
    mean_snapshot = np.mean(snapshots, axis=0)
    snapshots_centered = snapshots - mean_snapshot
    
    # SVD分解
    U, S, Vt = np.linalg.svd(snapshots_centered.T, full_matrices=False)
    
    # 计算能量捕获
    energy = np.cumsum(S**2) / np.sum(S**2)
    
    # 选择模态数（捕获99%能量）
    r = np.argmax(energy >= 0.99) + 1
    print(f"选择的模态数: {r} (捕获{energy[r-1]*100:.2f}%能量)")
    
    # 提取POD模态
    phi = U[:, :r]  # (nx*ny, r)
    
    # 计算模态系数
    coefficients = snapshots_centered @ phi  # (n_samples, r)
    
    # 构建代理模型（高斯过程）
    print("\n训练高斯过程代理模型...")
    from sklearn.gaussian_process import GaussianProcessRegressor
    from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C
    
    kernel = C(1.0, (1e-3, 1e3)) * RBF(1.0, (1e-2, 1e2))
    gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10)
    gp.fit(params, coefficients)
    
    print("代理模型训练完成！")
    
    # 测试降阶模型
    print("\n测试降阶模型...")
    n_test = 10
    test_params = np.random.uniform(0.5, 1.5, (n_test, 4))
    
    # 降阶模型预测
    pred_coefficients, sigma = gp.predict(test_params, return_std=True)
    pred_snapshots = pred_coefficients @ phi.T + mean_snapshot
    
    # 高保真求解对比
    hf_snapshots = []
    for i in range(n_test):
        T_left, T_right, T_top, T_bottom = test_params[i]
        
        T = np.zeros((ny, nx))
        T[:, 0] = T_left
        T[:, -1] = T_right
        T[0, :] = T_bottom
        T[-1, :] = T_top
        
        for _ in range(1000):
            T_old = T.copy()
            T[1:-1, 1:-1] = 0.25 * (T_old[1:-1, :-2] + T_old[1:-1, 2:] + 
                                     T_old[:-2, 1:-1] + T_old[2:, 1:-1])
            T[:, 0] = T_left
            T[:, -1] = T_right
            T[0, :] = T_bottom
            T[-1, :] = T_top
        
        hf_snapshots.append(T.flatten())
    
    hf_snapshots = np.array(hf_snapshots)
    
    # 计算误差
    errors = np.abs(pred_snapshots - hf_snapshots)
    mean_error = np.mean(errors, axis=1)
    max_error = np.max(errors, axis=1)
    
    print(f"平均误差: {np.mean(mean_error):.6f}")
    print(f"最大误差: {np.mean(max_error):.6f}")
    
    # 可视化
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    
    # 奇异值衰减
    axes[0, 0].semilogy(S[:20], 'bo-')
    axes[0, 0].axvline(x=r, color='r', linestyle='--', label=f'Selected r={r}')
    axes[0, 0].set_xlabel('Mode Index')
    axes[0, 0].set_ylabel('Singular Value')
    axes[0, 0].set_title('Singular Value Decay')
    axes[0, 0].legend()
    axes[0, 0].grid(True)
    
    # 能量捕获
    axes[0, 1].plot(energy[:20], 'g-')
    axes[0, 1].axvline(x=r, color='r', linestyle='--', label=f'99% Energy')
    axes[0, 1].axhline(y=0.99, color='k', linestyle=':', label='99% Threshold')
    axes[0, 1].set_xlabel('Number of Modes')
    axes[0, 1].set_ylabel('Cumulative Energy')
    axes[0, 1].set_title('Energy Captured by POD Modes')
    axes[0, 1].legend()
    axes[0, 1].grid(True)
    
    # 前4个POD模态
    for i in range(4):
        mode = phi[:, i].reshape(ny, nx)
        im = axes[0, 2].contourf(X, Y, mode, levels=20, cmap='RdBu_r')
        axes[0, 2].set_title(f'POD Mode {i+1}')
        axes[0, 2].set_xlabel('x')
        axes[0, 2].set_ylabel('y')
        plt.colorbar(im, ax=axes[0, 2])
        if i < 3:
            axes[0, 2] = plt.subplot(2, 3, 3)
    
    # 重新绘制模态
    fig2, axes2 = plt.subplots(2, 2, figsize=(10, 8))
    for i in range(4):
        ax = axes2[i//2, i%2]
        mode = phi[:, i].reshape(ny, nx)
        im = ax.contourf(X, Y, mode, levels=20, cmap='RdBu_r')
        ax.set_title(f'POD Mode {i+1}')
        ax.set_xlabel('x')
        ax.set_ylabel('y')
        plt.colorbar(im, ax=ax)
    plt.tight_layout()
    plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题043\\output\\case2_pod_modes.png', dpi=150)
    plt.close()
    
    # 预测vs真实对比
    test_idx = 0
    pred_field = pred_snapshots[test_idx].reshape(ny, nx)
    hf_field = hf_snapshots[test_idx].reshape(ny, nx)
    error_field = np.abs(pred_field - hf_field)
    
    im1 = axes[1, 0].contourf(X, Y, hf_field, levels=20, cmap='hot')
    axes[1, 0].set_title('High-Fidelity Solution')
    axes[1, 0].set_xlabel('x')
    axes[1, 0].set_ylabel('y')
    plt.colorbar(im1, ax=axes[1, 0])
    
    im2 = axes[1, 1].contourf(X, Y, pred_field, levels=20, cmap='hot')
    axes[1, 1].set_title('ROM Prediction')
    axes[1, 1].set_xlabel('x')
    axes[1, 1].set_ylabel('y')
    plt.colorbar(im2, ax=axes[1, 1])
    
    im3 = axes[1, 2].contourf(X, Y, error_field, levels=20, cmap='coolwarm')
    axes[1, 2].set_title(f'Error (Max: {np.max(error_field):.4f})')
    axes[1, 2].set_xlabel('x')
    axes[1, 2].set_ylabel('y')
    plt.colorbar(im3, ax=axes[1, 2])
    
    plt.tight_layout()
    plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题043\\output\\case2_pod_rom.png', dpi=150)
    plt.close()
    
    print("\n案例2完成！")
    return phi, coefficients, gp

# 运行案例2
if __name__ == "__main__":
    case2_pod_rom()

8.3 案例3：神经网络代理模型

def case3_nn_surrogate():
    """案例3：神经网络代理模型"""
    print("\n" + "="*60)
    print("案例3：神经网络代理模型")
    print("="*60)
    
    # 问题：悬臂梁端部位移预测
    # 输入：材料参数(E, rho)、几何参数(L, b, h)、载荷(P)
    # 输出：端部挠度、最大应力
    
    # 生成训练数据（解析解）
    print("\n生成训练数据...")
    n_samples = 1000
    
    # 参数范围
    E_range = [1e9, 3e11]  # 杨氏模量 (Pa)
    rho_range = [500, 10000]  # 密度 (kg/m^3)
    L_range = [0.5, 5.0]  # 长度 (m)
    b_range = [0.01, 0.5]  # 宽度 (m)
    h_range = [0.005, 0.3]  # 高度 (m)
    P_range = [100, 10000]  # 载荷 (N)
    
    # 拉丁超立方采样
    from pyDOE2 import lhs
    samples = lhs(6, samples=n_samples)
    
    E = E_range[0] + samples[:, 0] * (E_range[1] - E_range[0])
    rho = rho_range