主题077：数据驱动建模

目录

1. 引言
2. 机器学习基础
3. 神经网络在CFD中的应用
4. 物理信息神经网络(PINN)
5. 数据同化方法
6. 代理模型与降阶建模
7. Python仿真案例
8. 总结与展望

---

引言

1.1 数据驱动建模的兴起

传统的计算流体力学(CFD)方法基于第一性原理，通过求解Navier-Stokes方程等控制方程来模拟流动和传热现象。然而，随着计算需求的不断增长和数据获取能力的提升，一种全新的范式——数据驱动建模正在改变工程仿真的面貌。

什么是数据驱动建模？

数据驱动建模是指利用大量数据（实验数据、仿真数据或观测数据）来构建数学模型，而不是完全依赖物理定律的显式表达。这种方法的核心思想是：

• 从数据中学习系统的内在规律
• 将物理知识作为约束融入模型
• 实现快速预测和实时仿真

1.2 为什么需要数据驱动建模？

传统CFD面临的挑战：

1. 计算成本高：复杂的CFD模拟可能需要数天甚至数周的计算时间
2. 实时性不足：对于在线监测和实时控制，传统CFD无法满足时效要求
3. 参数化困难：多参数优化需要大量的重复计算
4. 不确定性量化：蒙特卡洛方法需要成千上万次模拟

数据驱动建模的优势：

1. 计算速度快：训练好的模型可以在毫秒级完成预测
2. 实时性强：适合在线监测和数字孪生应用
3. 参数化容易：可以快速评估不同参数组合
4. 融合多源数据：可以整合实验数据、仿真数据和观测数据

1.3 本教程内容概述

本教程将系统介绍数据驱动建模在对流换热仿真中的应用：

• 机器学习和深度学习基础
• 神经网络在CFD中的应用
• 物理信息神经网络(PINN)原理与实现
• 数据同化方法
• 代理模型与降阶建模技术
• 5个完整的Python仿真案例，包含GIF动画

---

机器学习基础

2.1 机器学习概述

机器学习是人工智能的一个分支，其核心是让计算机从数据中学习规律，而不是通过显式编程。根据学习方式的不同，机器学习可以分为：

监督学习：

• 从标记数据中学习输入-输出映射
• 示例：回归、分类
• CFD应用：从CFD数据学习流动特征

无监督学习：

• 从未标记数据中发现隐藏结构
• 示例：聚类、降维
• CFD应用：流动模态分析

强化学习：

• 通过与环境交互学习最优策略
• 示例：游戏AI、机器人控制
• CFD应用：流动控制优化

2.2 深度学习基础

深度学习是机器学习的一个子领域，使用多层神经网络来学习数据的层次化表示。

人工神经元：

人工神经元是神经网络的基本单元，其数学模型为：

$$y=f\left(\sum _{i=1}^n{w}_i{x}_i+b\right)$$

其中：

• $x_i$：输入特征
• $w_i$：权重
• $b$：偏置
• $f$：激活函数

常用激活函数：

激活函数	公式	特点
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	$\text{ReLU}(x)=\max (0,x)$	计算简单，缓解梯度消失
Leaky ReLU	$\max (\alpha x,x)$	解决ReLU死亡问题
Swish	$x\cdot \sigma (x)$	自门控，性能较好

前馈神经网络：

前馈神经网络（多层感知机）是最基本的神经网络结构：

$${\mathbf{z}}^{[l]}={\mathbf{W}}^{[l]}{\mathbf{a}}^{[l-1]}+{\mathbf{b}}^{[l]}$$

$${\mathbf{a}}^{[l]}=f({\mathbf{z}}^{[l]})$$

其中$l$表示层数。

2.3 训练神经网络

损失函数：

损失函数衡量模型预测与真实值之间的差距：

均方误差(MSE)：

$${\mathcal{L}}_{\text{MSE}}=\frac{1}{N}\sum _{i=1}^N(y_i-{\hat{y}}_i{)}^2$$

平均绝对误差(MAE)：

$${\mathcal{L}}_{\text{MAE}}=\frac{1}{N}\sum _{i=1}^N|y_i-{\hat{y}}_i|$$

交叉熵损失：

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

优化算法：

梯度下降：

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

随机梯度下降(SGD)：

使用小批量数据计算梯度，加快训练速度。

Adam优化器：

结合了动量和自适应学习率的优点：

$$m_t={\beta }_1{m}_{t-1}+(1-{\beta }_1)g_t$$

$$v_t={\beta }_2{v}_{t-1}+(1-{\beta }_2)g_t^2$$

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

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

正则化技术：

L2正则化（权重衰减）：

$${\mathcal{L}}_{\text{total}}=\mathcal{L}+\lambda \sum _i{w}_i^2$$

Dropout：

训练时随机丢弃一部分神经元，防止过拟合。

早停(Early Stopping)：

当验证集性能不再提升时停止训练。

2.4 深度学习框架

PyTorch：

• 动态计算图，易于调试
• 学术界广泛使用
• 适合研究和原型开发

TensorFlow：

• 静态计算图，部署方便
• 工业界应用广泛
• TensorBoard可视化强大

JAX：

• 自动微分、XLA编译
• 适合高性能科学计算
• 函数式编程风格

---

神经网络在CFD中的应用

3.1 流动特征提取

神经网络可以自动从流场数据中提取有意义的特征。

卷积神经网络(CNN)：

CNN特别适合处理网格化的流场数据：

$$(\mathbf{I}\ast \mathbf{K})(i,j)=\sum _m\sum _n\mathbf{I}(i+m,j+n)\mathbf{K}(m,n)$$

在CFD中的应用：

• 流场可视化分析
• 涡结构识别
• 流动模态分解

案例：涡识别：

使用CNN自动识别流场中的涡结构，替代传统的Q准则或λ2准则。

3.2 流场预测

全连接网络预测：

使用全连接网络建立参数到流场的映射：

$$\mathbf{u}(\mathbf{x};\theta )=\text{NN}(\mathbf{x};\mathbf{w})$$

卷积神经网络预测：

使用编码器-解码器结构：

输入参数 → 编码器 → 潜在空间 → 解码器 → 流场

U-Net架构：

U-Net在流场超分辨率和修正中表现出色：

• 编码器提取多尺度特征
• 跳跃连接保留细节信息
• 解码器重建高分辨率流场

3.3 湍流建模

雷诺应力建模：

使用神经网络替代传统的湍流模型：

$${\tau }_{ij}^{\text{RANS}}=\text{NN}(S_{ij},{\Omega }_{ij},\dots )$$

亚格子应力建模(LES)：

$${\tau }_{ij}^{\text{SGS}}=\text{NN}({\bar{S}}_{ij},\Delta ,\dots )$$

深度强化学习控制：

使用强化学习优化流动控制策略：

• 状态：流场特征
• 动作：控制输入（射流、壁面运动等）
• 奖励：减阻或增强换热

3.4 传热预测

换热系数预测：

建立几何参数和流动条件到换热系数的映射：

$$h=\text{NN}(Re,Pr,\text{geometry})$$

温度场预测：

预测非均匀边界条件下的温度分布：

$$T(\mathbf{x},t)=\text{NN}(\mathbf{x},t;\text{boundary conditions})$$

---

物理信息神经网络(PINN)

4.1 PINN的基本思想

物理信息神经网络(Physics-Informed Neural Networks, PINN)是将物理定律嵌入神经网络的一种新方法。

核心思想：

不仅用数据训练网络，还用物理方程约束网络：

$${\mathcal{L}}_{\text{total}}={\mathcal{L}}_{\text{data}}+{\lambda }_{\text{pde}}{\mathcal{L}}_{\text{pde}}+{\lambda }_{\text{bc}}{\mathcal{L}}_{\text{bc}}$$

其中：

• ${\mathcal{L}}_{\text{data}}$：数据拟合损失
• ${\mathcal{L}}_{\text{pde}}$：PDE残差损失
• ${\mathcal{L}}_{\text{bc}}$：边界条件损失

PINN的优势：

1. 数据效率高：即使数据稀疏也能获得合理结果
2. 物理一致性：解满足物理定律
3. 逆问题求解：可以反推未知参数
4. 无网格方法：不需要传统CFD的网格

4.2 PINN的数学框架

问题设定：

考虑一般形式的PDE：

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

带有边界条件：

$$\mathcal{B}[u](\mathbf{x},t)=0,\mathbf{x}\in \partial \Omega$$

和初始条件：

$$u(\mathbf{x},0)=u_0(\mathbf{x})$$

神经网络近似：

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

自动微分：

使用自动微分计算PDE残差：

$${\mathcal{L}}_{\text{pde}}=\frac{1}{N_f}\sum _{i=1}^{N_f}|\mathcal{F}[u_{\text{NN}}]({\mathbf{x}}_f^i,t_f^i){|}^2$$

PyTorch和TensorFlow都提供了自动微分功能：

# PyTorch示例
u_pred = model(x, t)
u_x = torch.autograd.grad(u_pred, x, grad_outputs=torch.ones_like(u_pred), create_graph=True)[0]
u_xx = torch.autograd.grad(u_x, x, grad_outputs=torch.ones_like(u_x), create_graph=True)[0]

4.3 对流扩散方程的PINN求解

控制方程：

$$\frac{\partial u}{\partial t}+\mathbf{v}\cdot \nabla u=\alpha {\nabla }^{2}u+f$$

PINN实现：

1. 定义神经网络$u_{\text{NN}}(x,y,t)$
2. 计算自动微分得到$u_t,u_x,u_y,u_{xx},u_{yy}$
3. 计算PDE残差：

$$r=u_t+v_x{u}_x+v_y{u}_y-\alpha (u_{xx}+u_{yy})-f$$

4. 损失函数：

$$\mathcal{L}=\frac{1}{N_u}\sum |u_{\text{pred}}-u_{\text{true}}{|}^2+\frac{1}{N_f}\sum |r{|}^2$$

4.4 Navier-Stokes方程的PINN求解

不可压NS方程：

连续性方程：

$$\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实现策略：

速度-压力联合求解：

定义网络输出：$[u,v,p]=\text{NN}(x,y,t)$

损失函数包含：

• 连续性方程残差
• x方向动量方程残差
• y方向动量方程残差
• 边界条件
• 观测数据（如有）

流函数-涡量形式：

对于二维流动，使用流函数$\psi$：

$$u=\frac{\partial \psi }{\partial y},v=-\frac{\partial \psi }{\partial x}$$

自动满足连续性方程，减少约束数量。

4.5 传热问题的PINN求解

能量方程：

$$\frac{\partial T}{\partial t}+\mathbf{u}\cdot \nabla T=\alpha {\nabla }^{2}T$$

与流动耦合：

1. 先求解NS方程得到速度场
2. 用速度场作为已知输入求解能量方程
3. 或者同时求解流动和温度场

Boussinesq近似：

对于自然对流，浮力项：

$$\mathbf{f}=-\beta (T-T_0)\mathbf{g}$$

4.6 PINN的改进与扩展

自适应权重：

不同损失项的重要性可能不同，使用自适应权重：

$${\lambda }_{\text{pde}}=\frac{|{\nabla }_{\theta }{\mathcal{L}}_{\text{data}}|}{|{\nabla }_{\theta }{\mathcal{L}}_{\text{pde}}|}$$

域分解PINN：

将计算域分解为多个子域，每个子域用一个网络：

• 减少单个网络的复杂度
• 提高训练效率
• 适合并行计算

时空自适应采样：

在残差大的区域增加采样点：

1. 在初始点集上训练
2. 评估PDE残差
3. 在残差大的区域添加新点
4. 重新训练

守恒律约束：

除了点-wise的PDE残差，还可以强制满足积分形式的守恒律。

---

数据同化方法

5.1 数据同化概述

数据同化是将观测数据与数值模型结合，以获得最优估计的方法。

基本问题：

给定：

• 数值模型预测：${\mathbf{x}}^f$（先验估计）
• 观测数据：${\mathbf{y}}^o$

求：最优估计${\mathbf{x}}^a$（后验估计）

应用背景：

• 天气预报
• 海洋环流模拟
• 实时流场重建
• 数字孪生

5.2 卡尔曼滤波

卡尔曼滤波(KF)：

对于线性系统：

状态方程：

$${\mathbf{x}}_k={\mathbf{F}}_k{\mathbf{x}}_{k-1}+{\mathbf{w}}_k$$

观测方程：

$${\mathbf{y}}_k={\mathbf{H}}_k{\mathbf{x}}_k+{\mathbf{v}}_k$$

其中${\mathbf{w}}_k$和${\mathbf{v}}_k$是过程噪声和观测噪声。

卡尔曼滤波步骤：

1. 预测：

$${\mathbf{x}}_k^f={\mathbf{F}}_k{\mathbf{x}}_{k-1}^a$$

$${\mathbf{P}}_k^f={\mathbf{F}}_k{\mathbf{P}}_{k-1}^a{\mathbf{F}}_k^T+{\mathbf{Q}}_k$$

2. 更新：

$${\mathbf{K}}_k={\mathbf{P}}_k^f{\mathbf{H}}_k^T({\mathbf{H}}_k{\mathbf{P}}_k^f{\mathbf{H}}_k^T+{\mathbf{R}}_k{)}^{-1}$$

$${\mathbf{x}}_k^a={\mathbf{x}}_k^f+{\mathbf{K}}_k({\mathbf{y}}_k-{\mathbf{H}}_k{\mathbf{x}}_k^f)$$

$${\mathbf{P}}_k^a=(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k){\mathbf{P}}_k^f$$

5.3 扩展卡尔曼滤波(EKF)

对于非线性系统：

状态方程：

$${\mathbf{x}}_k=f({\mathbf{x}}_{k-1})+{\mathbf{w}}_k$$

观测方程：

$${\mathbf{y}}_k=h({\mathbf{x}}_k)+{\mathbf{v}}_k$$

线性化：

使用雅可比矩阵：

$${\mathbf{F}}_k=\frac{\partial f}{\partial \mathbf{x}}{|}_{{\mathbf{x}}_{k-1}^a}$$

$${\mathbf{H}}_k=\frac{\partial h}{\partial \mathbf{x}}{|}_{{\mathbf{x}}_k^f}$$

5.4 集合卡尔曼滤波(EnKF)

基本思想：

用集合（ensemble）表示概率分布，避免计算协方差矩阵。

算法步骤：

1. 初始化：生成$N_e$个初始样本
2. 预报：每个样本独立积分
3. 更新：使用卡尔曼增益更新每个样本

优点：

• 适用于高维系统
• 不需要切线性模型
• 可以处理非高斯分布

缺点：

• 需要大量样本
• 样本间可能存在相关性

5.5 变分同化

三维变分同化(3D-Var)：

最小化代价函数：

$$J(\mathbf{x})=\frac{1}{2}(\mathbf{x}-{\mathbf{x}}^b{)}^T{\mathbf{B}}^{-1}(\mathbf{x}-{\mathbf{x}}^b)+\frac{1}{2}(\mathbf{y}-H(\mathbf{x}){)}^T{\mathbf{R}}^{-1}(\mathbf{y}-H(\mathbf{x}))$$

四维变分同化(4D-Var)：

考虑时间演化的变分同化：

$$J({\mathbf{x}}_0)=\frac{1}{2}({\mathbf{x}}_0-{\mathbf{x}}^b{)}^T{\mathbf{B}}^{-1}({\mathbf{x}}_0-{\mathbf{x}}^b)+\frac{1}{2}\sum _{i=0}^N({\mathbf{y}}_i-H({\mathbf{x}}_i){)}^T{\mathbf{R}}_i^{-1}({\mathbf{y}}_i-H({\mathbf{x}}_i))$$

其中${\mathbf{x}}_i=M_{0\to i}({\mathbf{x}}_0)$是模式积分结果。

5.6 粒子滤波

基本思想：

用一组带权重的粒子表示后验分布。

算法步骤：

1. 采样：从先验分布采样粒子
2. 预测：每个粒子按系统模型演化
3. 权重更新：根据观测计算权重
4. 重采样：根据权重重采样

优点：

• 可以处理任意分布
• 适用于强非线性系统

缺点：

• 维度灾难
• 需要大量粒子

---

代理模型与降阶建模

6.1 代理模型概述

代理模型（Surrogate Model）是复杂模型的廉价近似，用于快速预测。

应用场景：

• 设计优化（需要大量评估）
• 不确定性量化（蒙特卡洛模拟）
• 实时预测（在线应用）

常用代理模型：

方法	特点	适用场景
多项式响应面(RSM)	简单直观	低维问题
Kriging/Gaussian Process	提供不确定性估计	中等维度
径向基函数(RBF)	插值能力强	非线性问题
神经网络	表达能力强	高维复杂问题
支持向量机(SVM)	泛化能力强	分类和回归

6.2 高斯过程回归

高斯过程：

高斯过程是函数的分布，完全由均值函数和协方差函数决定：

$$f(\mathbf{x})\sim \mathcal{G}\mathcal{P}(m(\mathbf{x}),k(\mathbf{x},{\mathbf{x}}^{\prime }))$$

常用核函数：

平方指数核(RBF)：

$$k(\mathbf{x},{\mathbf{x}}^{\prime })={\sigma }^2\exp \left(-\frac{\Vert \mathbf{x}-{\mathbf{x}}^{\prime }{\Vert }^2}{2l^2}\right)$$

Matérn核：

$$k_{\nu }(\mathbf{x},{\mathbf{x}}^{\prime })={\sigma }^2\frac{2^{1-\nu }}{\Gamma (\nu )}{\left(\frac{\sqrt{2\nu }r}{l}\right)}^{\nu }{K}_{\nu }\left(\frac{\sqrt{2\nu }r}{l}\right)$$

预测：

给定训练数据$\mathbf{X},\mathbf{y}$，预测新点${\mathbf{x}}_{\ast }$：

$${\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 }$$

6.3 本征正交分解(POD)

POD基本原理：

寻找一组基函数，使得数据在这组基上的投影能量最大。

数学推导：

给定快照矩阵$\mathbf{S}=[{\mathbf{u}}_1,{\mathbf{u}}_2,...,{\mathbf{u}}_N]$，求解特征值问题：

$$\mathbf{R}\mathbf{\varphi }=\lambda \mathbf{\varphi }$$

其中$\mathbf{R}=\frac{1}{N}\mathbf{S}{\mathbf{S}}^T$是相关矩阵。

截断：

保留前$M$个模态，使得能量占比：

$$E=\frac{{\sum }_{i=1}^M{\lambda }_i}{{\sum }_{i=1}^N{\lambda }_i}>0.99$$

Galerkin投影：

将控制方程投影到POD模态上，得到降阶模型：

$$\frac{d\mathbf{a}}{dt}=\mathbf{Aa}+\mathbf{B}(\mathbf{a})$$

6.4 动态模态分解(DMD)

DMD基本原理：

从数据中提取动态特征，假设系统满足线性动力学：

$${\mathbf{x}}_{k+1}=\mathbf{A}{\mathbf{x}}_k$$

算法步骤：

1. 构建数据矩阵：$\mathbf{X}=[{\mathbf{x}}_1,...,{\mathbf{x}}_{N-1}]$，${\mathbf{X}}^{\prime }=[{\mathbf{x}}_2,...,{\mathbf{x}}_N]$
2. 计算近似线性算子：$\mathbf{A}={\mathbf{X}}^{\prime }{\mathbf{X}}^{+}$
3. 对$\mathbf{A}$进行特征值分解
4. 重构动力学

与POD的关系：

DMD可以看作是在POD基上的线性回归。

6.5 神经网络降阶模型

编码器-解码器结构：

高维流场 → 编码器 → 低维潜在空间 → 解码器 → 高维流场

自编码器降阶：

训练目标：

$$\underset{\theta }{\min }\Vert \mathbf{u}-\text{Decoder}(\text{Encoder}(\mathbf{u})){\Vert }^2$$

LSTM预测：

在潜在空间用LSTM预测时间演化：

$${\mathbf{z}}_{t+1}=\text{LSTM}({\mathbf{z}}_t)$$

6.6 物理约束的降阶模型

问题：纯数据驱动的降阶模型可能不满足物理约束。

解决方案：

1. 物理损失函数：

$$\mathcal{L}={\mathcal{L}}_{\text{reconstruction}}+\lambda {\mathcal{L}}_{\text{physics}}$$

2. 结构保持降阶：

保持原始系统的守恒性质。

3. 算子推断：

从数据推断低维算子，保持结构。

---

Python仿真案例

案例1：神经网络拟合非线性函数

本案例演示使用神经网络拟合复杂的非线性函数。

"""
案例1：神经网络拟合非线性函数
演示深度学习的基本原理
"""
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.animation as animation

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

print("案例1：神经网络拟合非线性函数")

# 生成训练数据：复杂的非线性函数
def true_function(x, y):
    """真实函数：复杂的非线性关系"""
    return np.sin(2*np.pi*x) * np.cos(2*np.pi*y) + 0.3*np.sin(4*np.pi*x*y)

# 生成数据
np.random.seed(42)
n_samples = 1000
x_train = np.random.uniform(0, 1, n_samples)
y_train = np.random.uniform(0, 1, n_samples)
z_train = true_function(x_train, y_train) + 0.1*np.random.randn(n_samples)

# 测试数据
x_test = np.linspace(0, 1, 50)
y_test = np.linspace(0, 1, 50)
X_test, Y_test = np.meshgrid(x_test, y_test)
Z_true = true_function(X_test, Y_test)

# 神经网络类
class NeuralNetwork:
    def __init__(self, input_size, hidden_size, output_size):
        # 初始化权重
        self.W1 = np.random.randn(input_size, hidden_size) * 0.1
        self.b1 = np.zeros((1, hidden_size))
        self.W2 = np.random.randn(hidden_size, hidden_size) * 0.1
        self.b2 = np.zeros((1, hidden_size))
        self.W3 = np.random.randn(hidden_size, output_size) * 0.1
        self.b3 = np.zeros((1, output_size))
        
        # 存储训练历史
        self.loss_history = []
        
    def relu(self, x):
        return np.maximum(0, x)
    
    def relu_derivative(self, x):
        return (x > 0).astype(float)
    
    def forward(self, X):
        # 前向传播
        self.z1 = np.dot(X, self.W1) + self.b1
        self.a1 = self.relu(self.z1)
        self.z2 = np.dot(self.a1, self.W2) + self.b2
        self.a2 = self.relu(self.z2)
        self.z3 = np.dot(self.a2, self.W3) + self.b3
        return self.z3
    
    def backward(self, X, y, learning_rate):
        m = X.shape[0]
        
        # 前向传播
        output = self.forward(X)
        
        # 计算损失
        loss = np.mean((output - y)**2)
        
        # 反向传播
        d_output = 2 * (output - y) / m
        
        dW3 = np.dot(self.a2.T, d_output)
        db3 = np.sum(d_output, axis=0, keepdims=True)
        
        d_a2 = np.dot(d_output, self.W3.T)
        d_z2 = d_a2 * self.relu_derivative(self.z2)
        dW2 = np.dot(self.a1.T, d_z2)
        db2 = np.sum(d_z2, axis=0, keepdims=True)
        
        d_a1 = np.dot(d_z2, self.W2.T)
        d_z1 = d_a1 * self.relu_derivative(self.z1)
        dW1 = np.dot(X.T, d_z1)
        db1 = np.sum(d_z1, axis=0, keepdims=True)
        
        # 更新权重
        self.W1 -= learning_rate * dW1
        self.b1 -= learning_rate * db1
        self.W2 -= learning_rate * dW2
        self.b2 -= learning_rate * db2
        self.W3 -= learning_rate * dW3
        self.b3 -= learning_rate * db3
        
        return loss
    
    def train(self, X, y, epochs, learning_rate):
        for epoch in range(epochs):
            loss = self.backward(X, y, learning_rate)
            self.loss_history.append(loss)
            if epoch % 1000 == 0:
                print(f"Epoch {epoch}, Loss: {loss:.6f}")

# 准备数据
X_train = np.column_stack([x_train, y_train])
y_train_reshaped = z_train.reshape(-1, 1)

# 创建和训练网络
nn = NeuralNetwork(input_size=2, hidden_size=64, output_size=1)
print("开始训练神经网络...")
nn.train(X_train, y_train_reshaped, epochs=10000, learning_rate=0.01)

# 预测
X_test_flat = np.column_stack([X_test.ravel(), Y_test.ravel()])
Z_pred = nn.forward(X_test_flat).reshape(X_test.shape)

# 可视化
fig = plt.figure(figsize=(16, 12))

# 真实函数
ax1 = fig.add_subplot(2, 3, 1, projection='3d')
ax1.plot_surface(X_test, Y_test, Z_true, cmap='viridis')
ax1.set_xlabel('X')
ax1.set_ylabel('Y')
ax1.set_zlabel('Z')
ax1.set_title('真实函数')

# 神经网络预测
ax2 = fig.add_subplot(2, 3, 2, projection='3d')
ax2.plot_surface(X_test, Y_test, Z_pred, cmap='viridis')
ax2.set_xlabel('X')
ax2.set_ylabel('Y')
ax2.set_zlabel('Z')
ax2.set_title('神经网络预测')

# 误差
ax3 = fig.add_subplot(2, 3, 3, projection='3d')
error = np.abs(Z_true - Z_pred)
ax3.plot_surface(X_test, Y_test, error, cmap='hot')
ax3.set_xlabel('X')
ax3.set_ylabel('Y')
ax3.set_zlabel('Error')
ax3.set_title('绝对误差')

# 训练数据分布
ax4 = fig.add_subplot(2, 3, 4)
scatter = ax4.scatter(x_train, y_train, c=z_train, cmap='viridis', s=10)
ax4.set_xlabel('X')
ax4.set_ylabel('Y')
ax4.set_title('训练数据分布')
plt.colorbar(scatter, ax=ax4)

# 损失曲线
ax5 = fig.add_subplot(2, 3, 5)
ax5.semilogy(nn.loss_history)
ax5.set_xlabel('Epoch')
ax5.set_ylabel('Loss (log scale)')
ax5.set_title('训练损失曲线')
ax5.grid(True)

# 预测vs真实值
ax6 = fig.add_subplot(2, 3, 6)
ax6.scatter(Z_true.ravel(), Z_pred.ravel(), alpha=0.5, s=1)
ax6.plot([-2, 2], [-2, 2], 'r--', label='Perfect prediction')
ax6.set_xlabel('True Value')
ax6.set_ylabel('Predicted Value')
ax6.set_title('预测准确性')
ax6.legend()
ax6.grid(True)

plt.tight_layout()
plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题077_数据驱动建模\\仿真结果\\case1_neural_network_fitting.png', dpi=150, bbox_inches='tight')
plt.close()

print("案例1完成！")

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

本案例演示使用物理信息神经网络求解热传导方程。

"""
案例2：PINN求解一维热传导方程
物理信息神经网络入门示例
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

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

print("案例2：PINN求解一维热传导方程")

# 物理参数
alpha = 0.1  # 热扩散系数
L = 1.0      # 域长度
T_final = 1.0  # 最终时间

# 神经网络类（简化版）
class PINN:
    def __init__(self, layers):
        self.layers = layers
        self.weights = []
        self.biases = []
        
        # 初始化权重
        for i in range(len(layers) - 1):
            w = np.random.randn(layers[i], layers[i+1]) * np.sqrt(2.0 / layers[i])
            b = np.zeros((1, layers[i+1]))
            self.weights.append(w)
            self.biases.append(b)
        
        self.loss_history = {'total': [], 'ic': [], 'bc': [], 'pde': []}
    
    def tanh(self, x):
        return np.tanh(x)
    
    def tanh_derivative(self, x):
        return 1 - np.tanh(x)**2
    
    def forward(self, x, t):
        # 合并输入
        X = np.hstack([x, t])
        
        # 前向传播
        self.activations = [X]
        for i in range(len(self.weights) - 1):
            z = np.dot(self.activations[-1], self.weights[i]) + self.biases[i]
            a = self.tanh(z)
            self.activations.append(a)
        
        # 输出层（线性激活）
        z = np.dot(self.activations[-1], self.weights[-1]) + self.biases[-1]
        self.activations.append(z)
        
        return z
    
    def compute_derivatives(self, x, t):
        """计算自动微分"""
        eps = 1e-5
        
        # u
        u = self.forward(x, t)
        
        # u_t
        u_t = (self.forward(x, t + eps) - self.forward(x, t - eps)) / (2 * eps)
        
        # u_x
        u_x = (self.forward(x + eps, t) - self.forward(x - eps, t)) / (2 * eps)
        
        # u_xx
        u_xx = (self.forward(x + eps, t) - 2*u + self.forward(x - eps, t)) / (eps**2)
        
        return u, u_t, u_x, u_xx
    
    def compute_loss(self, x_ic, t_ic, u_ic, x_bc, t_bc, u_bc, x_f, t_f):
        """计算总损失"""
        # 初始条件损失
        u_ic_pred = self.forward(x_ic, t_ic)
        loss_ic = np.mean((u_ic_pred - u_ic)**2)
        
        # 边界条件损失
        u_bc_pred = self.forward(x_bc, t_bc)
        loss_bc = np.mean((u_bc_pred - u_bc)**2)
        
        # PDE残差损失
        u_f, u_t, u_x, u_xx = self.compute_derivatives(x_f, t_f)
        f = u_t - alpha * u_xx
        loss_pde = np.mean(f**2)
        
        # 总损失
        loss_total = loss_ic + loss_bc + loss_pde
        
        return loss_total, loss_ic, loss_bc, loss_pde
    
    def train_step(self, x_ic, t_ic, u_ic, x_bc, t_bc, u_bc, x_f, t_f, lr):
        """单步训练（数值梯度）"""
        eps = 1e-5
        
        # 计算当前损失
        loss_total, loss_ic, loss_bc, loss_pde = self.compute_loss(
            x_ic, t_ic, u_ic, x_bc, t_bc, u_bc, x_f, t_f
        )
        
        # 更新权重（简化版梯度下降）
        for i in range(len(self.weights)):
            # 权重梯度
            for j in range(self.weights[i].shape[0]):
                for k in range(self.weights[i].shape[1]):
                    self.weights[i][j, k] += eps
                    loss_plus, _, _, _ = self.compute_loss(
                        x_ic, t_ic, u_ic, x_bc, t_bc, u_bc, x_f, t_f
                    )
                    grad = (loss_plus - loss_total) / eps
                    self.weights[i][j, k] -= eps + lr * grad
            
            # 偏置梯度
            for j in range(self.biases[i].shape[1]):
                self.biases[i][0, j] += eps
                loss_plus, _, _, _ = self.compute_loss(
                    x_ic, t_ic, u_ic, x_bc, t_bc, u_bc, x_f, t_f
                )
                grad = (loss_plus - loss_total) / eps
                self.biases[i][0, j] -= eps + lr * grad
        
        return loss_total, loss_ic, loss_bc, loss_pde
    
    def train(self, x_ic, t_ic, u_ic, x_bc, t_bc, u_bc, x_f, t_f, epochs, lr):
        """训练网络"""
        print("开始训练PINN...")
        for epoch in range(epochs):
            loss_total, loss_ic, loss_bc, loss_pde = self.train_step(
                x_ic, t_ic, u_ic, x_bc, t_bc, u_bc, x_f, t_f, lr
            )
            
            self.loss_history['total'].append(loss_total)
            self.loss_history['ic'].append(loss_ic)
            self.loss_history['bc'].append(loss_bc)
            self.loss_history['pde'].append(loss_pde)
            
            if epoch % 1000 == 0:
                print(f"Epoch {epoch}: Total={loss_total:.6f}, IC={loss_ic:.6f}, BC={loss_bc:.6f}, PDE={loss_pde:.6f}")

# 生成训练数据
# 初始条件：u(x, 0) = sin(pi*x)
x_ic = np.random.rand(100, 1) * L
t_ic = np.zeros((100, 1))
u_ic = np.sin(np.pi * x_ic)

# 边界条件：u(0, t) = u(1, t) = 0
x_bc = np.vstack([np.zeros((50, 1)), np.ones((50, 1)) * L])
t_bc = np.random.rand(100, 1) * T_final
u_bc = np.zeros((100, 1))

# 内部配点
x_f = np.random.rand(1000, 1) * L
t_f = np.random.rand(1000, 1) * T_final

# 创建和训练PINN
pinn = PINN(layers=[2, 32, 32, 32, 1])
pinn.train(x_ic, t_ic, u_ic, x_bc, t_bc, u_bc, x_f, t_f, epochs=5000, lr=0.001)

# 解析解
def analytical_solution(x, t, alpha, n_terms=20):
    """一维热传导方程的解析解"""
    u = np.zeros_like(x)
    for n in range(1, n_terms + 1):
        u += np.sin(n * np.pi * x) * np.exp(-(n * np.pi)**2 * alpha * t)
    u *= 2 / np.pi
    return u

# 生成预测网格
x_pred = np.linspace(0, L, 100)
t_pred = np.linspace(0, T_final, 50)
X_pred, T_pred = np.meshgrid(x_pred, t_pred)

# PINN预测
U_pred = np.zeros_like(X_pred)
for i in range(X_pred.shape[0]):
    for j in range(X_pred.shape[1]):
        U_pred[i, j] = pinn.forward(
            np.array([[X_pred[i, j]]]), 
            np.array([[T_pred[i, j]]])
        )[0, 0]

# 解析解
U_analytical = analytical_solution(X_pred, T_pred, alpha)

# 可视化
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

# PINN解
ax = axes[0, 0]
im = ax.contourf(X_pred, T_pred, U_pred, levels=20, cmap='hot')
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_title('PINN预测')
plt.colorbar(im, ax=ax)

# 解析解
ax = axes[0, 1]
im = ax.contourf(X_pred, T_pred, U_analytical, levels=20, cmap='hot')
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_title('解析解')
plt.colorbar(im, ax=ax)

# 误差
ax = axes[0, 2]
error = np.abs(U_pred - U_analytical)
im = ax.contourf(X_pred, T_pred, error, levels=20, cmap='viridis')
ax.set_xlabel('x')
ax.set_ylabel('t')
ax.set_title('绝对误差')
plt.colorbar(im, ax=ax)

# 损失曲线
ax = axes[1, 0]
ax.semilogy(pinn.loss_history['total'], label='Total', linewidth=2)
ax.semilogy(pinn.loss_history['ic'], label='IC', alpha=0.7)
ax.semilogy(pinn.loss_history['bc'], label='BC', alpha=0.7)
ax.semilogy(pinn.loss_history['pde'], label='PDE', alpha=0.7)
ax.set_xlabel('Epoch')
ax.set_ylabel('Loss (log scale)')
ax.set_title('训练损失曲线')
ax.legend()
ax.grid(True)

# 不同时刻的温度分布
ax = axes[1, 1]
times = [0, 0.2, 0.5, 1.0]
for t_val in times:
    idx = np.argmin(np.abs(t_pred - t_val))
    ax.plot(x_pred, U_pred[idx, :], label=f't={t_val}')
ax.set_xlabel('x')
ax.set_ylabel('u')
ax.set_title('不同时刻的温度分布')
ax.legend()
ax.grid(True)

# 空间-时间演化
ax = axes[1, 2]
ax.plot(x_pred, U_analytical[-1, :], 'b-', label='Analytical', linewidth=2)
ax.plot(x_pred, U_pred[-1, :], 'r--', label='PINN', linewidth=2)
ax.set_xlabel('x')
ax.set_ylabel('u')
ax.set_title(f't={T_final}时刻对比')
ax.legend()
ax.grid(True)

plt.tight_layout()
plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题077_数据驱动建模\\仿真结果\\case2_pinn_heat_equation.png', dpi=150, bbox_inches='tight')
plt.close()

print("案例2完成！")

案例3：高斯过程回归构建代理模型

本案例演示使用高斯过程回归构建计算昂贵的CFD模型的代理模型。

"""
案例3：高斯过程回归构建代理模型
使用GP替代昂贵的CFD计算
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist

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

print("案例3：高斯过程回归构建代理模型")

# 高斯过程回归类
class GaussianProcess:
    def __init__(self, kernel='rbf', noise=1e-5):
        self.kernel_type = kernel
        self.noise = noise
        self.X_train = None
        self.y_train = None
        self.K = None
        self.L = None
        self.alpha = None
        
    def rbf_kernel(self, X1, X2, length_scale=1.0, sigma=1.0):
        """RBF核函数"""
        sqdist = cdist(X1, X2, 'sqeuclidean')
        return sigma**2 * np.exp(-0.5 * sqdist / length_scale**2)
    
    def fit(self, X, y, length_scale=1.0, sigma=1.0):
        """训练GP模型"""
        self.X_train = X
        self.y_train = y
        self.length_scale = length_scale
        self.sigma = sigma
        
        # 计算核矩阵
        self.K = self.rbf_kernel(X, X, length_scale, sigma)
        self.K += self.noise * np.eye(len(X))
        
        # Cholesky分解
        self.L = np.linalg.cholesky(self.K)
        
        # 求解
        self.alpha = np.linalg.solve(self.L.T, np.linalg.solve(self.L, y))
        
    def predict(self, X_test):
        """预测"""
        # 计算测试点与训练点的核
        K_s = self.rbf_kernel(self.X_train, X_test, self.length_scale, self.sigma)
        
        # 预测均值
        mu = np.dot(K_s.T, self.alpha)
        
        # 预测方差
        K_ss = self.rbf_kernel(X_test, X_test, self.length_scale, self.sigma)
        v = np.linalg.solve(self.L, K_s)
        var = np.diag(K_ss) - np.sum(v**2, axis=0)
        std = np.sqrt(np.maximum(var, 0))
        
        return mu, std

# 模拟昂贵的CFD函数：努塞尔数关联式
def expensive_cfd_simulation(Re, Pr):
    """
    模拟昂贵的CFD计算：管内湍流换热
    使用Gnielinski关联式作为"真实"CFD
    """
    # Gnielinski关联式
    f = (0.79 * np.log(Re) - 1.64)**(-2)
    Nu = (f/8) * (Re - 1000) * Pr / (1 + 12.7 * np.sqrt(f/8) * (Pr**(2/3) - 1))
    # 添加一些"CFD噪声"
    Nu += 0.01 * Nu * np.random.randn(*Nu.shape)
    return Nu

# 生成训练数据（少量昂贵的CFD点）
np.random.seed(42)
n_train = 20
Re_train = np.random.uniform(10000, 100000, n_train)
Pr_train = np.random.uniform(0.7, 10, n_train)
X_train = np.column_stack([np.log10(Re_train), np.log10(Pr_train)])
y_train = expensive_cfd_simulation(Re_train, Pr_train)

# 创建GP模型
gp = GaussianProcess()
print("训练GP模型...")
gp.fit(X_train, y_train, length_scale=0.5, sigma=50)

# 生成测试数据（密集网格）
Re_test = np.linspace(10000, 100000, 100)
Pr_test = np.linspace(0.7, 10, 100)
Re_grid, Pr_grid = np.meshgrid(Re_test, Pr_test)
X_test = np.column_stack([np.log10(Re_grid.ravel()), np.log10(Pr_grid.ravel())])

# GP预测
print("进行GP预测...")
mu, std = gp.predict(X_test)
mu = mu.reshape(Re_grid.shape)
std = std.reshape(Re_grid.shape)

# "真实"CFD（用于对比）
Nu_true = expensive_cfd_simulation(Re_grid, Pr_grid)

# 可视化
fig, axes = plt.subplots(2, 3, figsize=(16, 10))

# GP预测均值
ax = axes[0, 0]
im = ax.contourf(Re_grid, Pr_grid, mu, levels=20, cmap='viridis')
ax.scatter(Re_train, Pr_train, c='red', s=50, marker='x', label='Training points')
ax.set_xlabel('Re')
ax.set_ylabel('Pr')
ax.set_title('GP预测均值')
ax.set_xscale('log')
plt.colorbar(im, ax=ax)
ax.legend()

# 真实CFD
ax = axes[0, 1]
im = ax.contourf(Re_grid, Pr_grid, Nu_true, levels=20, cmap='viridis')
ax.set_xlabel('Re')
ax.set_ylabel('Pr')
ax.set_title('真实CFD结果')
ax.set_xscale('log')
plt.colorbar(im, ax=ax)

# 预测标准差（不确定性）
ax = axes[0, 2]
im = ax.contourf(Re_grid, Pr_grid, std, levels=20, cmap='hot')
ax.scatter(Re_train, Pr_train, c='blue', s=50, marker='x')
ax.set_xlabel('Re')
ax.set_ylabel('Pr')
ax.set_title('预测不确定性')
ax.set_xscale('log')
plt.colorbar(im, ax=ax)

# 误差分布
ax = axes[1, 0]
error = np.abs(mu - Nu_true)
im = ax.contourf(Re_grid, Pr_grid, error, levels=20, cmap='Reds')
ax.set_xlabel('Re')
ax.set_ylabel('Pr')
ax.set_title('绝对误差')
ax.set_xscale('log')
plt.colorbar(im, ax=ax)

# 1D切片对比（固定Pr=5）
ax = axes[1, 1]
Pr_fixed = 5.0
idx = np.argmin(np.abs(Pr_test - Pr_fixed))
ax.plot(Re_test, Nu_true[idx, :], 'b-', label='True CFD', linewidth=2)
ax.plot(Re_test, mu[idx, :], 'r--', label='GP Prediction', linewidth=2)
ax.fill_between(Re_test, 
                mu[idx, :] - 2*std[idx, :], 
                mu[idx, :] + 2*std[idx, :], 
                alpha=0.3, color='red', label='95% CI')
ax.scatter(Re_train[Pr_train > 4], y_train[Pr_train > 4], 
           c='green', s=50, marker='o', label='Training data')
ax.set_xlabel('Re')
ax.set_ylabel('Nu')
ax.set_title(f'Pr={Pr_fixed}时的预测对比')
ax.set_xscale('log')
ax.legend()
ax.grid(True)

# 误差统计
ax = axes[1, 2]
errors = (mu - Nu_true).ravel()
ax.hist(errors, bins=30, edgecolor='black', alpha=0.7)
ax.axvline(0, color='r', linestyle='--', linewidth=2)
ax.set_xlabel('Prediction Error')
ax.set_ylabel('Frequency')
ax.set_title(f'误差分布 (RMSE={np.sqrt(np.mean(errors**2)):.2f})')
ax.grid(True)

plt.tight_layout()
plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题077_数据驱动建模\\仿真结果\\case3_gaussian_process.png', dpi=150, bbox_inches='tight')
plt.close()

print("案例3完成！")

案例4：本征正交分解(POD)降阶模型

本案例演示使用POD从CFD数据中提取主要模态并构建降阶模型。

"""
案例4：本征正交分解(POD)降阶模型
从CFD数据中提取主要流动模态
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import eigh

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

print("案例4：本征正交分解(POD)降阶模型")

# 生成模拟CFD数据：圆柱绕流的不同时间步
# 简化模型：使用解析函数模拟卡门涡街

def generate_flow_field(x, y, t, Re=100):
    """生成模拟流场数据"""
    nx, ny = len(x), len(y)
    X, Y = np.meshgrid(x, y, indexing='ij')
    
    # 圆柱位置
    cx, cy = 2.0, 2.0
    
    # 基础流动 + 涡脱落
    u = np.ones_like(X) * 0.5
    v = np.zeros_like(Y)
    
    # 添加涡结构
    for i, ti in enumerate(t):
        St = 0.2  # 斯特劳哈尔数
        omega = 2 * np.pi * St * 0.5 / 0.5  # 涡脱落频率
        
        # 涡对
        for j in range(3):
            x_vortex = cx + (j+1) * 1.5 + 0.5 * ti
            y_vortex_p = cy + 0.3 * np.sin(omega * ti + j * np.pi)
            y_vortex_n = cy - 0.3 * np.sin(omega * ti + j * np.pi)
            
            r_p = np.sqrt((X - x_vortex)**2 + (Y - y_vortex_p)**2)
            r_n = np.sqrt((X - x_vortex)**2 + (Y - y_vortex_n)**2)
            
            strength = 0.3 * np.exp(-0.1 * (j+1))
            u += strength * (Y - y_vortex_p) / (r_p + 0.1)**2 * np.exp(-r_p/2)
            u -= strength * (Y - y_vortex_n) / (r_n + 0.1)**2 * np.exp(-r_n/2)
            v -= strength * (X - x_vortex) / (r_p + 0.1)**2 * np.exp(-r_p/2)
            v += strength * (X - x_vortex) / (r_n + 0.1)**2 * np.exp(-r_n/2)
    
    return u, v

# 生成快照数据
nx, ny = 60, 40
x = np.linspace(0, 10, nx)
y = np.linspace(0, 4, ny)
n_snapshots = 100
t_snapshots = np.linspace(0, 20, n_snapshots)

print(f"生成{n_snapshots}个流场快照...")
snapshots = []
for i, t in enumerate(t_snapshots):
    u, v = generate_flow_field(x, y, [t])
    # 合并u和v作为状态向量
    state = np.concatenate([u.flatten(), v.flatten()])
    snapshots.append(state)

snapshots = np.array(snapshots).T  # 形状: (2*nx*ny, n_snapshots)

# POD分析
print("执行POD分析...")

# 计算均值并去均值
mean_flow = np.mean(snapshots, axis=1, keepdims=True)
snapshots_centered = snapshots - mean_flow

# 计算相关矩阵
C = np.dot(snapshots_centered.T, snapshots_centered) / n_snapshots

# 特征值分解
eigenvalues, eigenvectors = eigh(C)

# 按特征值降序排序
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]

# 计算POD模态
pod_modes = np.dot(snapshots_centered, eigenvectors)

# 归一化
for i in range(pod_modes.shape[1]):
    norm = np.sqrt(np.sum(pod_modes[:, i]**2))
    if norm > 0:
        pod_modes[:, i] /= norm

# 计算能量占比
energy = eigenvalues / np.sum(eigenvalues)
cumulative_energy = np.cumsum(energy)

# 可视化
fig, axes = plt.subplots(2, 3, figsize=(16, 10))

# 能量分布
ax = axes[0, 0]
ax.semilogy(range(1, 21), energy[:20], 'bo-', markersize=6)
ax.set_xlabel('Mode Number')
ax.set_ylabel('Energy')
ax.set_title('POD模态能量分布')
ax.grid(True)

# 累积能量
ax = axes[0, 1]
ax.plot(range(1, 21), cumulative_energy[:20], 'r-', linewidth=2)
ax.axhline(0.95, color='g', linestyle='--', label='95% Energy')
ax.axhline(0.99, color='b', linestyle='--', label='99% Energy')
ax.set_xlabel('Number of Modes')
ax.set_ylabel('Cumulative Energy')
ax.set_title('累积能量占比')
ax.legend()
ax.grid(True)

# 前4个POD模态
for i in range(4):
    ax = axes[0 if i < 2 else 1, 2 if i % 2 == 0 else 0 if i == 1 else 1]
    mode = pod_modes[:, i]
    u_mode = mode[:nx*ny].reshape(nx, ny)
    
    im = ax.contourf(x, y, u_mode.T, levels=20, cmap='RdBu_r')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_title(f'POD Mode {i+1} (Energy: {energy[i]*100:.1f}%)')
    plt.colorbar(im, ax=ax)

# 重构流场（使用前10个模态）
n_modes = 10
ax = axes[1, 2]

# 选择最后一个快照进行重构
snapshot_idx = -1
coefficients = np.dot(snapshots_centered[:, snapshot_idx], pod_modes[:, :n_modes])
reconstructed = mean_flow[:, 0] + np.dot(pod_modes[:, :n_modes], coefficients)

u_original = snapshots[:nx*ny, snapshot_idx].reshape(nx, ny)
u_reconstructed = reconstructed[:nx*ny].reshape(nx, ny)

# 显示重构误差
error = np.abs(u_original - u_reconstructed)
im = ax.contourf(x, y, error.T, levels=20, cmap='hot')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title(f'重构误差 ({n_modes} modes)')
plt.colorbar(im, ax=ax)

plt.tight_layout()
plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题077_数据驱动建模\\仿真结果\\case4_pod_analysis.png', dpi=150, bbox_inches='tight')
plt.close()

# 绘制重构对比图
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

ax = axes[0]
im = ax.contourf(x, y, u_original.T, levels=20, cmap='RdBu_r')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Original Flow')
plt.colorbar(im, ax=ax)

ax = axes[1]
im = ax.contourf(x, y, u_reconstructed.T, levels=20, cmap='RdBu_r')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title(f'Reconstructed ({n_modes} modes)')
plt.colorbar(im, ax=ax)

ax = axes[2]
im = ax.contourf(x, y, error.T, levels=20, cmap='hot')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Reconstruction Error')
plt.colorbar(im, ax=ax)

plt.tight_layout()
plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题077_数据驱动建模\\仿真结果\\case4_pod_reconstruction.png', dpi=150, bbox_inches='tight')
plt.close()

print("案例4完成！")

案例5：数据同化融合实验与仿真数据

本案例演示使用卡尔曼滤波融合实验数据和CFD仿真结果。

"""
案例5：数据同化融合实验与仿真数据
使用卡尔曼滤波融合多源数据
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

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

print("案例5：数据同化融合实验与仿真数据")

# 一维热传导问题的数据同化
# 真实系统（未知）
def true_system(x, t, alpha=0.1):
    """真实温度场"""
    return np.sin(np.pi * x) * np.exp(-alpha * (np.pi)**2 * t)

# CFD模型（有偏差）
def cfd_model(x, t, alpha_cfd=0.08):
    """CFD模型（参数有偏差）"""
    return np.sin(np.pi * x) * np.exp(-alpha_cfd * (np.pi)**2 * t)

# 生成网格
nx = 50
x = np.linspace(0, 1, nx)
dx = x[1] - x[0]

# 时间设置
dt = 0.01
nt = 100
t = np.arange(nt) * dt

# 初始化
T_true = np.zeros((nt, nx))
T_cfd = np.zeros((nt, nx))
T_da = np.zeros((nt, nx))  # 数据同化结果

# 初始条件
T_true[0, :] = np.sin(np.pi * x)
T_cfd[0, :] = np.sin(np.pi * x)
T_da[0, :] = np.sin(np.pi * x)

# 观测设置（稀疏观测）
n_obs = 5
obs_indices = np.linspace(5, nx-5, n_obs, dtype=int)
obs_noise = 0.05

# 卡尔曼滤波参数
# 状态转移矩阵（简化的一维热传导离散）
alpha_da = 0.09  # 数据同化使用的参数
r = alpha_da * dt / dx**2
F = np.eye(nx)
for i in range(1, nx-1):
    F[i, i-1] = r
    F[i, i] = 1 - 2*r
    F[i, i+1] = r

# 观测矩阵（只观测特定位置）
H = np.zeros((n_obs, nx))
for i, idx in enumerate(obs_indices):
    H[i, idx] = 1

# 协方差矩阵
Q = 0.001 * np.eye(nx)  # 过程噪声
R = obs_noise**2 * np.eye(n_obs)  # 观测噪声
P = 0.01 * np.eye(nx)  # 初始误差协方差

print("运行数据同化...")

# 存储历史
obs_history = []
rmse_cfd = []
rmse_da = []

for n in range(1, nt):
    # 真实系统演化
    T_true[n, :] = true_system(x, t[n])
    
    # CFD模型预测（有偏差）
    T_cfd[n, :] = cfd_model(x, t[n])
    
    # 数据同化：预测步骤
    T_pred = F @ T_da[n-1, :]
    P_pred = F @ P @ F.T + Q
    
    # 生成观测（稀疏 + 噪声）
    if n % 5 == 0:  # 每5步观测一次
        obs = T_true[n, obs_indices] + obs_noise * np.random.randn(n_obs)
        obs_history.append((n, obs))
        
        # 更新步骤
        K = P_pred @ H.T @ np.linalg.inv(H @ P_pred @ H.T + R)
        T_da[n, :] = T_pred + K @ (obs - H @ T_pred)
        P = (np.eye(nx) - K @ H) @ P_pred
    else:
        T_da[n, :] = T_pred
        P = P_pred
    
    # 计算RMSE
    rmse_cfd.append(np.sqrt(np.mean((T_cfd[n, :] - T_true[n, :])**2)))
    rmse_da.append(np.sqrt(np.mean((T_da[n, :] - T_true[n, :])**2)))

# 生成GIF动画
print("生成GIF动画...")
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

# 初始化图像
lines = []
for ax in axes:
    line1, = ax.plot([], [], 'b-', label='True', linewidth=2)
    line2, = ax.plot([], [], 'r--', label='CFD', linewidth=2)
    line3, = ax.plot([], [], 'g:', label='Data Assimilation', linewidth=2)
    lines.extend([line1, line2, line3])

# 设置坐标轴
for ax in axes:
    ax.set_xlim(0, 1)
    ax.set_ylim(-0.2, 1.2)
    ax.set_xlabel('x')
    ax.set_ylabel('Temperature')
    ax.legend()
    ax.grid(True)

axes[0].set_title('Temperature Distribution')
axes[1].set_title('Error (CFD vs True)')
axes[2].set_title('Error (DA vs True)')

def init():
    for line in lines:
        line.set_data([], [])
    return lines

def update(frame):
    n = frame
    
    # 温度分布
    lines[0].set_data(x, T_true[n, :])
    lines[1].set_data(x, T_cfd[n, :])
    lines[2].set_data(x, T_da[n, :])
    
    # CFD误差
    error_cfd = T_cfd[n, :] - T_true[n, :]
    lines[3].set_data([], [])
    lines[4].set_data(x, error_cfd)
    lines[5].set_data([], [])
    
    # DA误差
    error_da = T_da[n, :] - T_true[n, :]
    lines[6].set_data([], [])
    lines[7].set_data([], [])
    lines[8].set_data(x, error_da)
    
    axes[0].set_title(f'Temperature Distribution (t={t[n]:.2f})')
    
    return lines

# 创建动画
anim = FuncAnimation(fig, update, frames=range(0, nt, 2), init_func=init, blit=True)
anim.save('d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题077_数据驱动建模\\仿真结果\\case5_data_assimilation.gif', 
          writer='pillow', fps=10, dpi=100)
plt.close()

# 绘制最终结果
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 最后时刻对比
ax = axes[0, 0]
ax.plot(x, T_true[-1, :], 'b-', label='True', linewidth=2)
ax.plot(x, T_cfd[-1, :], 'r--', label='CFD', linewidth=2)
ax.plot(x, T_da[-1, :], 'g:', label='Data Assimilation', linewidth=2)
ax.scatter(x[obs_indices], T_true[-1, obs_indices], c='orange', s=50, marker='o', label='Observations')
ax.set_xlabel('x')
ax.set_ylabel('Temperature')
ax.set_title(f'Final Time (t={t[-1]:.2f})')
ax.legend()
ax.grid(True)

# RMSE演化
ax = axes[0, 1]
ax.plot(t[1:], rmse_cfd, 'r-', label='CFD RMSE', linewidth=2)
ax.plot(t[1:], rmse_da, 'g-', label='DA RMSE', linewidth=2)
ax.set_xlabel('Time')
ax.set_ylabel('RMSE')
ax.set_title('RMSE Evolution')
ax.legend()
ax.grid(True)

# 误差分布
ax = axes[1, 0]
error_cfd_final = T_cfd[-1, :] - T_true[-1, :]
error_da_final = T_da[-1, :] - T_true[-1, :]
ax.plot(x, error_cfd_final, 'r-', label='CFD Error', linewidth=2)
ax.plot(x, error_da_final, 'g-', label='DA Error', linewidth=2)
ax.axhline(0, color='k', linestyle='--', alpha=0.5)
ax.set_xlabel('x')
ax.set_ylabel('Error')
ax.set_title('Final Error Distribution')
ax.legend()
ax.grid(True)

# 误差统计
ax = axes[1, 1]
ax.bar(['CFD', 'Data Assimilation'], 
       [np.mean(rmse_cfd), np.mean(rmse_da)], 
       color=['red', 'green'], alpha=0.7)
ax.set_ylabel('Mean RMSE')
ax.set_title('Average RMSE Comparison')
ax.grid(True, axis='y')

plt.tight_layout()
plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题077_数据驱动建模\\仿真结果\\case5_data_assimilation_summary.png', dpi=150, bbox_inches='tight')
plt.close()

print("案例5完成！")

---