主题017：基于机器学习的结构优化

一、引言

1.1 背景与动机

结构优化设计作为工程领域的重要研究方向，其核心目标是在满足各种约束条件的前提下，寻找最优的材料分布或几何形状，以实现结构性能的最大化。传统的结构优化方法，如尺寸优化、形状优化和拓扑优化，虽然在理论上已经相当成熟，但在实际工程应用中仍面临诸多挑战。

首先，传统优化方法通常需要进行大量的有限元分析（FEA）计算。以拓扑优化为例，每一次迭代都需要重新组装刚度矩阵、求解线性方程组，计算灵敏度信息。对于复杂的三维结构或精细的网格划分，单次有限元分析可能需要数分钟甚至数小时。当优化问题涉及数百甚至数千次迭代时，总的计算时间可能达到数天或数周，这对于需要快速响应的设计迭代过程来说是不可接受的。

其次，传统优化方法对于设计空间的探索能力有限。许多优化算法，如梯度下降法、序列线性规划法等，都是基于局部搜索策略，容易陷入局部最优解。虽然全局优化算法如遗传算法、模拟退火算法等可以一定程度上缓解这一问题，但这些方法的计算成本更高，且收敛速度较慢。

近年来，机器学习（Machine Learning, ML）技术的快速发展为解决上述问题提供了新的思路。机器学习算法，特别是深度学习（Deep Learning）方法，具有强大的非线性映射能力和特征提取能力，可以从大量的数据中学习复杂的模式，并建立高效的预测模型。将机器学习技术应用于结构优化领域，可以在以下几个方面带来显著的改进：

1. 加速优化过程：通过训练代理模型（Surrogate Model）来替代昂贵的有限元分析，可以大幅减少优化过程中的计算时间。

2. 提高设计质量：机器学习算法可以从大量的历史优化案例中学习设计规律，帮助工程师发现更优的设计方案。

3. 实现实时优化：训练好的机器学习模型可以在毫秒级别内给出优化结果，使得实时交互式设计成为可能。

4. 处理高维设计空间：深度学习模型可以有效处理高维设计空间，对于传统方法难以解决的复杂优化问题具有优势。
   
   

1.2 学习目标

通过本主题的学习，读者将能够：

• 理解机器学习在结构优化中的基本原理和应用场景
• 掌握代理模型的构建方法和训练技巧
• 学会使用神经网络进行拓扑结构预测
• 了解基于机器学习的快速优化框架
• 能够独立开发简单的机器学习辅助优化程序

1.3 应用场景

基于机器学习的结构优化技术在以下领域具有广泛的应用前景：

• 航空航天：飞机机翼、发动机叶片等复杂结构的轻量化设计
• 汽车工业：车身框架、底盘结构的多目标优化
• 建筑工程：大跨度桥梁、高层建筑的抗震优化设计
• 生物医学：人工关节、骨科植入物的个性化设计
• 增材制造：考虑制造工艺约束的拓扑优化

二、核心理论

2.1 机器学习基础

2.1.1 监督学习与无监督学习

机器学习算法主要分为监督学习（Supervised Learning）和无监督学习（Unsupervised Learning）两大类。

监督学习是指从标注数据中学习映射函数的过程。在结构优化中，监督学习的典型应用包括：

• 回归问题：预测结构的性能指标（如应力、位移、固有频率等）
• 分类问题：判断结构是否满足设计要求
• 生成问题：根据边界条件和载荷生成最优拓扑结构

监督学习的数学表述为：给定训练数据集 $D=\{(x_i,y_i){\}}_{i=1}^n$，其中 $x_i$ 是输入特征，$y_i$ 是对应的标签，目标是学习一个映射函数 $f:X\to Y$，使得预测误差最小化：

$$\underset{f}{\min }\sum _{i=1}^{n}L(f(x_i),y_i)+\lambda R(f)$$

其中，$L$ 是损失函数，$R(f)$ 是正则化项，$\lambda$ 是正则化系数。

无监督学习是指从未标注数据中发现隐藏模式的过程。在结构优化中，无监督学习可以用于：

• 聚类分析：将相似的结构设计方案分组
• 降维：将高维设计空间映射到低维空间，便于可视化分析
• 特征学习：自动提取结构特征，用于后续优化

2.1.2 神经网络基础

神经网络是机器学习中最强大的模型之一，其核心思想是模拟生物神经系统的信息处理方式。一个典型的神经网络由输入层、隐藏层和输出层组成，每层包含多个神经元（节点）。

前向传播：给定输入 $x$，神经网络的输出通过逐层计算得到：

$$z^{[l]}=W^{[l]}{a}^{[l-1]}+b^{[l]}$$
$$a^{[l]}=g(z^{[l]})$$

其中，$W^{[l]}$ 和 $b^{[l]}$ 分别是第 $l$ 层的权重矩阵和偏置向量，$g(\cdot )$ 是激活函数，常用的激活函数包括：

• Sigmoid函数：$g(z)=\frac{1}{1+e^{-z}}$
• ReLU函数：$g(z)=\max (0,z)$
• Tanh函数：$g(z)=\frac{e^z-e^{-z}}{e^z+e^{-z}}$

反向传播：通过链式法则计算损失函数对网络参数的梯度，然后使用梯度下降法更新参数：

$$W^{[l]}=W^{[l]}-\alpha \frac{\partial L}{\partial W^{[l]}}$$
$$b^{[l]}=b^{[l]}-\alpha \frac{\partial L}{\partial b^{[l]}}$$

其中，$\alpha$ 是学习率。

2.1.3 卷积神经网络（CNN）

卷积神经网络是一种专门用于处理网格结构数据（如图像）的神经网络。在结构优化中，拓扑结构可以自然地表示为二维或三维的密度场，因此CNN特别适合用于拓扑结构的预测和生成。

CNN的核心组件包括：

卷积层：通过卷积核（滤波器）提取局部特征：

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

其中，$I$ 是输入特征图，$K$ 是卷积核。

池化层：降低特征图的空间维度，减少计算量：

• 最大池化：$y={\max }_{i,j\in R}{x}_{i,j}$
• 平均池化：$y=\frac{1}{|R|}{\sum }_{i,j\in R}{x}_{i,j}$

全连接层：将卷积层提取的特征映射到输出空间。

2.2 代理模型方法

2.2.1 代理模型的概念

代理模型（Surrogate Model），也称为元模型（Metamodel）或响应面模型（Response Surface Model），是一种用于近似复杂计算模型的数学模型。在结构优化中，有限元分析通常计算成本高昂，代理模型可以在保证一定精度的前提下，大幅降低计算时间。

代理模型的构建过程通常包括以下步骤：

1. 试验设计（Design of Experiments, DOE）：在设计空间中选取有代表性的样本点
2. 样本计算：在样本点上运行高保真模型（如有限元分析）获得响应值
3. 模型训练：使用样本数据训练代理模型
4. 模型验证：使用独立的测试数据验证代理模型的精度

2.2.2 常用的代理模型

多项式响应面（Polynomial Response Surface, PRS）：

$$\hat{y}(x)={\beta }_0+\sum _{i=1}^n{\beta }_i{x}_i+\sum _{i=1}^n\sum _{j=i}^n{\beta }_{ij}{x}_i{x}_j+\cdots$$

多项式响应面模型简单直观，但对于高度非线性的问题精度有限。

克里金模型（Kriging Model）：

克里金模型是一种基于高斯过程的插值方法，不仅可以给出预测值，还可以给出预测的不确定性：

$$\hat{y}(x)=\mu +Z(x)$$

其中，$\mu$ 是全局趋势，$Z(x)$ 是均值为零的高斯随机过程。

径向基函数（Radial Basis Function, RBF）：

$$\hat{y}(x)=\sum _{i=1}^n{w}_i\varphi (\Vert x-x_i\Vert )$$

其中，$\varphi (\cdot )$ 是径向基函数，常用的包括高斯函数、多二次函数等。

神经网络代理模型：

如前所述，神经网络可以作为强大的代理模型，特别是对于高维、高度非线性的问题。

2.2.3 代理模型的精度评估

评估代理模型精度的常用指标包括：

均方根误差（Root Mean Square Error, RMSE）：

$$RMSE=\sqrt{\frac{1}{n}\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2}$$

决定系数（$R^2$）：

$$R^2=1-\frac{{\sum }_{i=1}^n(y_i-{\hat{y}}_i{)}^2}{{\sum }_{i=1}^n(y_i-\bar{y}{)}^2}$$

最大绝对误差（Maximum Absolute Error, MAE）：

$$MAE=\underset{i}{\max }|y_i-{\hat{y}}_i|$$

2.3 基于机器学习的拓扑优化

2.3.1 数据驱动的方法

数据驱动的拓扑优化方法通过训练机器学习模型来学习从边界条件、载荷到最优拓扑的映射关系。这种方法的核心思想是：一旦模型训练完成，对于新的设计问题，可以直接通过前向传播获得优化结果，无需进行迭代优化。

数据驱动方法的优点：

• 推理速度快，可实现实时优化
• 不需要迭代，避免了收敛性问题
• 可以学习到隐含的设计规律

数据驱动方法的缺点：

• 需要大量的训练数据
• 泛化能力受限于训练数据的覆盖范围
• 对于训练数据之外的问题可能表现不佳

2.3.2 物理信息神经网络（Physics-Informed Neural Networks, PINNs）

PINNs是一种将物理定律嵌入神经网络损失函数的方法，可以在缺乏大量标注数据的情况下训练神经网络。对于结构优化问题，可以将平衡方程、本构关系等物理约束作为软约束加入损失函数：

$$L=L_{data}+{\lambda }_1{L}_{physics}+{\lambda }_2{L}_{boundary}$$

其中，$L_{data}$ 是数据损失，$L_{physics}$ 是物理约束损失，$L_{boundary}$ 是边界条件损失。

2.3.3 强化学习在结构优化中的应用

强化学习（Reinforcement Learning, RL）是一种通过与环境交互来学习最优策略的机器学习方法。在结构优化中，可以将优化过程建模为一个马尔可夫决策过程（MDP）：

• 状态（State）：当前的材料分布
• 动作（Action）：材料密度的调整
• 奖励（Reward）：结构性能的改善
• 策略（Policy）：决定在给定状态下采取什么动作

强化学习的优势在于可以学习到全局最优策略，且不需要大量的标注数据。但强化学习的训练通常需要大量的交互，计算成本较高。

三、案例实战

3.1 案例一：基于神经网络的拓扑结构预测

本案例展示如何使用神经网络学习从载荷条件到最优拓扑的映射关系。

3.1.1 问题描述

考虑一个二维悬臂梁结构，左端固定，右端承受不同方向的集中载荷。目标是训练一个神经网络，能够根据载荷的大小和方向，预测最优的材料分布。

3.1.2 数据生成

首先，需要生成训练数据。对于每个样本，随机生成载荷条件，然后运行传统的SIMP拓扑优化算法获得最优拓扑。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation

# 设置matplotlib后端
plt.switch_backend('Agg')

class TopologyOptimizer:
    """
    SIMP拓扑优化器
    
    用于生成训练数据
    """
    
    def __init__(self, nelx=60, nely=20, volfrac=0.4, penal=3.0, rmin=1.5):
        self.nelx = nelx
        self.nely = nely
        self.volfrac = volfrac
        self.penal = penal
        self.rmin = rmin
        self.ndof = 2 * (nelx + 1) * (nely + 1)
        
        # 初始化密度场
        self.x = np.ones((nely, nelx)) * volfrac
        
        # 准备有限元
        self._prepare_fem()
        self._prepare_filter()
    
    def _prepare_fem(self):
        """准备有限元分析"""
        E0, nu = 1.0, 0.3
        
        # 单元刚度矩阵
        k = np.array([
            [1/2-nu/6, 1/8+nu/8, -1/4-nu/12, -1/8+3*nu/8, -1/4+nu/12, -1/8-nu/8, nu/6, 1/8-3*nu/8],
            [1/8+nu/8, 1/2-nu/6, -1/8+3*nu/8, nu/6, -1/8-nu/8, -1/4+nu/12, 1/8-3*nu/8, -1/4-nu/12],
            [-1/4-nu/12, -1/8+3*nu/8, 1/2-nu/6, -1/8-nu/8, nu/6, 1/8-3*nu/8, -1/4+nu/12, 1/8+nu/8],
            [-1/8+3*nu/8, nu/6, -1/8-nu/8, 1/2-nu/6, 1/8-3*nu/8, -1/4-nu/12, 1/8+nu/8, -1/4+nu/12],
            [-1/4+nu/12, -1/8-nu/8, nu/6, 1/8-3*nu/8, 1/2-nu/6, 1/8+nu/8, -1/4-nu/12, -1/8+3*nu/8],
            [-1/8-nu/8, -1/4+nu/12, 1/8-3*nu/8, -1/4-nu/12, 1/8+nu/8, 1/2-nu/6, -1/8+3*nu/8, nu/6],
            [nu/6, 1/8-3*nu/8, -1/4+nu/12, 1/8+nu/8, -1/4-nu/12, -1/8+3*nu/8, 1/2-nu/6, -1/8-nu/8],
            [1/8-3*nu/8, -1/4-nu/12, 1/8+nu/8, -1/4+nu/12, -1/8+3*nu/8, nu/6, -1/8-nu/8, 1/2-nu/6]
        ])
        
        self.KE = k / (1 - nu**2)
        
        # 节点编号
        nodenrs = np.arange((self.nelx+1)*(self.nely+1)).reshape(self.nely+1, self.nelx+1)
        edofVec = 2*nodenrs[0:-1,0:-1].reshape(self.nelx*self.nely,1)+1
        edofMat = np.tile(edofVec,(1,8)) + np.tile(
            np.array([0, 1, 2*self.nely+2, 2*self.nely+3, 2*self.nely, 2*self.nely+1, -2, -1]),
            (self.nelx*self.nely,1)
        )
        self.edofMat = edofMat
        
        # 组装索引
        iK = np.tile(np.reshape(edofMat,(-1,1)),(1,8))
        jK = np.tile(np.reshape(edofMat,(1,-1)),(8,1))
        self.iK = iK.reshape(-1)
        self.jK = jK.reshape(-1)
    
    def _prepare_filter(self):
        """准备密度过滤"""
        nfilter = int(self.nelx * self.nely * ((2*(np.ceil(self.rmin)-1)+1)**2))
        iH = np.zeros(nfilter)
        jH = np.zeros(nfilter)
        sH = np.zeros(nfilter)
        cc = 0
        
        for i in range(self.nelx):
            for j in range(self.nely):
                row = i * self.nely + j
                kk1 = int(np.maximum(i - (np.ceil(self.rmin) - 1), 0))
                kk2 = int(np.minimum(i + np.ceil(self.rmin), self.nelx))
                ll1 = int(np.maximum(j - (np.ceil(self.rmin) - 1), 0))
                ll2 = int(np.minimum(j + np.ceil(self.rmin), self.nely))
                
                for k in range(kk1, kk2):
                    for l in range(ll1, ll2):
                        col = k * self.nely + l
                        fac = self.rmin - np.sqrt(((i-k)*(i-k)+(j-l)*(j-l)))
                        iH[cc] = row
                        jH[cc] = col
                        sH[cc] = np.maximum(0.0, fac)
                        cc = cc + 1
        
        H = np.sparse.coo_matrix((sH,(iH,jH)),shape=(self.nelx*self.nely,self.nelx*self.nely)).tocsc()
        Hs = H.sum(1)
        self.H = H
        self.Hs = Hs
    
    def optimize(self, load_x, load_y, n_iter=100):
        """
        运行拓扑优化
        
        参数:
            load_x, load_y: 载荷方向分量
            n_iter: 迭代次数
        
        返回:
            x: 最优密度场
            history: 优化历史
        """
        # 应用载荷
        F = np.zeros(self.ndof)
        F[2*(self.nelx+1)*(self.nely+1)-1] = load_y  # 右下角y方向
        F[2*(self.nelx+1)*(self.nely+1)-2] = load_x  # 右下角x方向
        
        # 固定左边界
        fixeddofs = np.arange(2*(self.nely+1))
        freedofs = np.setdiff1d(np.arange(self.ndof), fixeddofs)
        
        history = {'compliance': [], 'volume': []}
        
        for iter in range(n_iter):
            # 组装刚度矩阵
            Emin, E0 = 1e-9, 1.0
            sK = np.tile(self.KE.flatten(), self.nelx*self.nely) * \
                 np.reshape((Emin + (self.x.flatten()**self.penal)*(E0-Emin)),(-1,1)).repeat(8,axis=1).flatten()
            K = np.sparse.coo_matrix((sK,(self.iK,self.jK)),shape=(self.ndof,self.ndof)).tocsc()
            K = (K + K.T) / 2
            
            # 求解
            U = np.zeros(self.ndof)
            try:
                U[freedofs] = np.linalg.solve(K[freedofs,:][:,freedofs].toarray(), F[freedofs])
            except:
                U[freedofs] = 0
            
            # 计算柔度
            ce = np.zeros(self.nelx*self.nely)
            for el in range(self.nelx*self.nely):
                ue = U[self.edofMat[el,:]]
                ce[el] = ue @ self.KE @ ue
            
            ce = ce.reshape(self.nely, self.nelx)
            c = np.sum((Emin + self.x**self.penal*(E0-Emin))*ce)
            
            # 灵敏度分析
            dc = -self.penal*(E0-Emin)*self.x**(self.penal-1)*ce
            dv = np.ones((self.nely, self.nelx))
            
            # 过滤灵敏度
            dc = self._filter_sensitivity(dc)
            
            # OC更新
            self.x = self._oc_update(dc, dv)
            
            history['compliance'].append(c)
            history['volume'].append(np.mean(self.x))
        
        return self.x.copy(), history
    
    def _filter_sensitivity(self, dc):
        """过滤灵敏度"""
        dc_flat = dc.flatten('F')
        dc_filtered = (self.H @ (dc_flat * self.x.flatten('F')) / self.Hs).reshape(self.nely, self.nelx, order='F')
        return dc_filtered
    
    def _oc_update(self, dc, dv):
        """OC优化准则更新"""
        l1, l2 = 0, 1e9
        move = 0.2
        xnew = np.zeros_like(self.x)
        
        while (l2 - l1) / (l1 + l2) > 1e-3:
            lmid = 0.5 * (l2 + l1)
            xnew = np.maximum(0.0, np.maximum(self.x - move, 
                          np.minimum(1.0, np.minimum(self.x + move, 
                          self.x * np.sqrt(-dc / dv / lmid)))))
            
            if np.sum(xnew) > self.volfrac * self.nelx * self.nely:
                l1 = lmid
            else:
                l2 = lmid
        
        return xnew


class SimpleNeuralNetwork:
    """
    简化神经网络（纯NumPy实现）
    
    用于拓扑预测的代理模型
    """
    
    def __init__(self, input_size, hidden_sizes, output_size, learning_rate=0.001):
        """
        初始化神经网络
        
        参数:
            input_size: 输入层大小
            hidden_sizes: 隐藏层大小列表
            output_size: 输出层大小
            learning_rate: 学习率
        """
        self.layers = []
        self.learning_rate = learning_rate
        self.loss_history = []
        
        # 构建网络层
        sizes = [input_size] + hidden_sizes + [output_size]
        for i in range(len(sizes) - 1):
            # Xavier初始化
            W = np.random.randn(sizes[i], sizes[i+1]) * np.sqrt(2.0 / sizes[i])
            b = np.zeros((1, sizes[i+1]))
            self.layers.append({'W': W, 'b': b})
    
    def relu(self, x):
        """ReLU激活函数"""
        return np.maximum(0, x)
    
    def relu_derivative(self, x):
        """ReLU导数"""
        return (x > 0).astype(float)
    
    def sigmoid(self, x):
        """Sigmoid激活函数"""
        return 1 / (1 + np.exp(-np.clip(x, -500, 500)))
    
    def forward(self, X):
        """
        前向传播
        
        参数:
            X: 输入数据 (n_samples, input_size)
        
        返回:
            output: 网络输出
            cache: 缓存中间结果用于反向传播
        """
        cache = []
        a = X
        
        for i, layer in enumerate(self.layers):
            z = a @ layer['W'] + layer['b']
            cache.append({'z': z, 'a': a})
            
            if i < len(self.layers) - 1:  # 隐藏层使用ReLU
                a = self.relu(z)
            else:  # 输出层使用Sigmoid
                a = self.sigmoid(z)
        
        return a, cache
    
    def compute_loss(self, y_true, y_pred):
        """计算均方误差损失"""
        return np.mean((y_true - y_pred)**2)
    
    def backward(self, X, y_true, y_pred, cache):
        """
        反向传播
        
        参数:
            X: 输入数据
            y_true: 真实标签
            y_pred: 预测输出
            cache: 前向传播的缓存
        """
        n_samples = X.shape[0]
        
        # 输出层梯度
        dz = (y_pred - y_true) * y_pred * (1 - y_pred)  # Sigmoid导数
        
        for i in range(len(self.layers) - 1, -1, -1):
            a_prev = cache[i]['a']
            
            # 计算梯度
            dW = (a_prev.T @ dz) / n_samples
            db = np.mean(dz, axis=0, keepdims=True)
            
            # 更新参数
            self.layers[i]['W'] -= self.learning_rate * dW
            self.layers[i]['b'] -= self.learning_rate * db
            
            # 传播到前一层
            if i > 0:
                da = dz @ self.layers[i]['W'].T
                dz = da * self.relu_derivative(cache[i-1]['z'])
    
    def train(self, X, y, epochs=100, batch_size=32, validation_data=None, verbose=1):
        """
        训练网络
        
        参数:
            X: 训练输入
            y: 训练输出
            epochs: 训练轮数
            batch_size: 批次大小
            validation_data: 验证数据 (X_val, y_val)
            verbose: 是否打印进度
        """
        n_samples = X.shape[0]
        n_batches = int(np.ceil(n_samples / batch_size))
        
        for epoch in range(epochs):
            # 打乱数据
            indices = np.random.permutation(n_samples)
            X_shuffled = X[indices]
            y_shuffled = y[indices]
            
            epoch_loss = 0
            
            for i in range(n_batches):
                start_idx = i * batch_size
                end_idx = min((i + 1) * batch_size, n_samples)
                
                X_batch = X_shuffled[start_idx:end_idx]
                y_batch = y_shuffled[start_idx:end_idx]
                
                # 前向传播
                output, cache = self.forward(X_batch)
                
                # 计算损失
                batch_loss = self.compute_loss(y_batch, output)
                epoch_loss += batch_loss * (end_idx - start_idx)
                
                # 反向传播
                self.backward(X_batch, y_batch, output, cache)
            
            epoch_loss /= n_samples
            self.loss_history.append(epoch_loss)
            
            if verbose and (epoch + 1) % 10 == 0:
                msg = f"Epoch {epoch+1}/{epochs}, Loss: {epoch_loss:.6f}"
                if validation_data:
                    X_val, y_val = validation_data
                    val_pred, _ = self.forward(X_val)
                    val_loss = self.compute_loss(y_val, val_pred)
                    msg += f", Val Loss: {val_loss:.6f}"
                print(msg)
    
    def predict(self, X):
        """预测"""
        output, _ = self.forward(X)
        return output


def generate_training_data(n_samples=100, nelx=30, nely=10):
    """
    生成训练数据
    
    参数:
        n_samples: 样本数量
        nelx, nely: 网格尺寸
    
    返回:
        X: 输入特征 (载荷条件)
        Y: 输出标签 (最优拓扑)
    """
    print(f"生成{n_samples}个训练样本...")
    
    optimizer = TopologyOptimizer(nelx=nelx, nely=nely, volfrac=0.4)
    
    X = []
    Y = []
    
    for i in range(n_samples):
        # 随机生成载荷
        angle = np.random.uniform(-np.pi/2, np.pi/2)
        magnitude = np.random.uniform(0.5, 2.0)
        load_x = magnitude * np.cos(angle)
        load_y = -abs(magnitude * np.sin(angle))  # 向下
        
        # 运行优化
        x_opt, _ = optimizer.optimize(load_x, load_y, n_iter=50)
        
        # 存储数据
        X.append([load_x, load_y])
        Y.append(x_opt.flatten())
        
        if (i + 1) % 10 == 0:
            print(f"  已生成 {i+1}/{n_samples} 个样本")
    
    return np.array(X), np.array(Y)


def main():
    """主程序"""
    print("="*60)
    print("基于机器学习的结构优化 - 神经网络拓扑预测")
    print("="*60)
    
    # 设置随机种子
    np.random.seed(42)
    
    # 生成训练数据
    n_samples = 200
    nelx, nely = 30, 10
    X, Y = generate_training_data(n_samples, nelx, nely)
    
    # 划分训练集和测试集
    split_idx = int(0.8 * n_samples)
    X_train, X_test = X[:split_idx], X[split_idx:]
    Y_train, Y_test = Y[:split_idx], Y[split_idx:]
    
    print(f"\n训练集大小: {len(X_train)}")
    print(f"测试集大小: {len(X_test)}")
    
    # 创建神经网络
    input_size = 2  # 载荷x, y分量
    hidden_sizes = [128, 256, 128]
    output_size = nelx * nely
    
    print(f"\n网络结构: {input_size} -> {hidden_sizes} -> {output_size}")
    
    nn = SimpleNeuralNetwork(input_size, hidden_sizes, output_size, learning_rate=0.01)
    
    # 训练网络
    print("\n开始训练神经网络...")
    nn.train(X_train, Y_train, epochs=200, batch_size=16, 
             validation_data=(X_test, Y_test), verbose=1)
    
    # 评估模型
    print("\n评估模型性能...")
    Y_pred = nn.predict(X_test)
    mse = np.mean((Y_test - Y_pred)**2)
    print(f"测试集MSE: {mse:.6f}")
    
    # 可视化结果
    print("\n生成可视化结果...")
    fig, axes = plt.subplots(3, 3, figsize=(12, 10))
    
    for i in range(3):
        idx = i * 5
        
        # 真实拓扑
        ax = axes[i, 0]
        true_topo = Y_test[idx].reshape(nely, nelx)
        ax.imshow(1 - true_topo, cmap='gray', vmin=0, vmax=1)
        ax.set_title(f'真实拓扑 (载荷: [{X_test[idx][0]:.2f}, {X_test[idx][1]:.2f}])')
        ax.axis('off')
        
        # 预测拓扑
        ax = axes[i, 1]
        pred_topo = Y_pred[idx].reshape(nely, nelx)
        ax.imshow(1 - pred_topo, cmap='gray', vmin=0, vmax=1)
        ax.set_title('神经网络预测')
        ax.axis('off')
        
        # 误差
        ax = axes[i, 2]
        error = np.abs(true_topo - pred_topo)
        im = ax.imshow(error, cmap='hot', vmin=0, vmax=1)
        ax.set_title(f'绝对误差 (MAE: {np.mean(error):.4f})')
        ax.axis('off')
        plt.colorbar(im, ax=ax, fraction=0.046)
    
    plt.tight_layout()
    plt.savefig('ml_comparison.png', dpi=150, bbox_inches='tight')
    print("对比图已保存为 ml_comparison.png")
    
    # 绘制训练历史
    fig, ax = plt.subplots(figsize=(10, 5))
    ax.plot(nn.loss_history, linewidth=2)
    ax.set_xlabel('Epoch')
    ax.set_ylabel('Loss (MSE)')
    ax.set_title('神经网络训练历史')
    ax.grid(True, alpha=0.3)
    ax.set_yscale('log')
    plt.tight_layout()
    plt.savefig('training_history.png', dpi=150, bbox_inches='tight')
    print("训练历史图已保存为 training_history.png")
    
    print("\n" + "="*60)
    print("程序执行完成！")
    print("="*60)


if __name__ == "__main__":
    main()

3.1.3 代码解析

数据生成部分：

def generate_training_data(n_samples=100, nelx=30, nely=10):
    """生成训练数据"""
    optimizer = TopologyOptimizer(nelx=nelx, nely=nely, volfrac=0.4)
    
    X = []
    Y = []
    
    for i in range(n_samples):
        # 随机生成载荷
        angle = np.random.uniform(-np.pi/2, np.pi/2)
        magnitude = np.random.uniform(0.5, 2.0)
        load_x = magnitude * np.cos(angle)
        load_y = -abs(magnitude * np.sin(angle))
        
        # 运行优化
        x_opt, _ = optimizer.optimize(load_x, load_y, n_iter=50)
        
        X.append([load_x, load_y])
        Y.append(x_opt.flatten())

这段代码首先创建一个拓扑优化器，然后循环生成多个训练样本。对于每个样本，随机生成载荷的大小和方向，运行SIMP优化算法获得最优拓扑，并将载荷条件作为输入特征、最优拓扑作为输出标签存储。

神经网络架构：

class SimpleNeuralNetwork:
    def __init__(self, input_size, hidden_sizes, output_size, learning_rate=0.001):
        self.layers = []
        sizes = [input_size] + hidden_sizes + [output_size]
        for i in range(len(sizes) - 1):
            W = np.random.randn(sizes[i], sizes[i+1]) * np.sqrt(2.0 / sizes[i])
            b = np.zeros((1, sizes[i+1]))
            self.layers.append({'W': W, 'b': b})

这里使用纯NumPy实现了一个简单的全连接神经网络。网络采用Xavier初始化策略，即权重初始化为均值为0、方差为 $2/n_{in}$ 的正态分布，这有助于加速训练收敛。

前向传播：

def forward(self, X):
    cache = []
    a = X
    
    for i, layer in enumerate(self.layers):
        z = a @ layer['W'] + layer['b']
        cache.append({'z': z, 'a': a})
        
        if i < len(self.layers) - 1:
            a = self.relu(z)
        else:
            a = self.sigmoid(z)
    
    return a, cache

前向传播依次通过每一层，计算线性变换 $z=aW+b$，然后应用激活函数。隐藏层使用ReLU激活函数，输出层使用Sigmoid激活函数将输出限制在[0, 1]范围内，这与密度场的物理意义相符。

反向传播：

def backward(self, X, y_true, y_pred, cache):
    n_samples = X.shape[0]
    dz = (y_pred - y_true) * y_pred * (1 - y_pred)
    
    for i in range(len(self.layers) - 1, -1, -1):
        a_prev = cache[i]['a']
        dW = (a_prev.T @ dz) / n_samples
        db = np.mean(dz, axis=0, keepdims=True)
        
        self.layers[i]['W'] -= self.learning_rate * dW
        self.layers[i]['b'] -= self.learning_rate * db
        
        if i > 0:
            da = dz @ self.layers[i]['W'].T
            dz = da * self.relu_derivative(cache[i-1]['z'])

反向传播使用链式法则计算损失函数对各层参数的梯度。对于Sigmoid输出层，梯度为 $(y_{pred}-y_{true})\cdot y_{pred}\cdot (1-y_{pred})$。然后使用梯度下降法更新参数。

3.1.4 运行结果预期

运行上述代码后，将生成以下结果：

1. 训练过程输出：显示每10个epoch的训练损失和验证损失，观察损失值逐渐下降
2. ml_comparison.png：对比图显示真实拓扑、神经网络预测拓扑和预测误差
3. training_history.png：训练损失随epoch变化的曲线

预期神经网络能够在测试集上达到较好的预测精度，预测结果与真实拓扑在视觉上相似，主要结构特征能够被正确捕捉。

3.2 案例二：代理模型加速优化

本案例展示如何使用代理模型替代昂贵的有限元分析，加速结构优化过程。

3.2.1 问题描述

考虑一个桁架结构的尺寸优化问题。目标是最小化结构重量，同时满足应力和位移约束。传统的优化方法需要在每次迭代中进行有限元分析，计算成本较高。本案例将使用神经网络代理模型来预测结构的响应，从而加速优化过程。

3.2.2 完整代码实现

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

plt.switch_backend('Agg')


class TrussFEM:
    """桁架有限元分析器"""
    
    def __init__(self, nodes, elements, E=2.1e11, rho=7850):
        """
        初始化桁架结构
        
        参数:
            nodes: 节点坐标 (n_nodes, 2)
            elements: 单元连接 (n_elements, 2)
            E: 弹性模量
            rho: 材料密度
        """
        self.nodes = nodes
        self.elements = elements
        self.E = E
        self.rho = rho
        self.n_nodes = len(nodes)
        self.n_elements = len(elements)
        self.ndof = 2 * self.n_nodes
    
    def analyze(self, areas, loads, fixed_dofs):
        """
        进行有限元分析
        
        参数:
            areas: 单元截面积
            loads: 节点载荷
            fixed_dofs: 固定自由度
        
        返回:
            displacement: 节点位移
            stress: 单元应力
            weight: 结构重量
        """
        # 组装刚度矩阵
        K = np.zeros((self.ndof, self.ndof))
        
        for i, (n1, n2) in enumerate(self.elements):
            x1, y1 = self.nodes[n1]
            x2, y2 = self.nodes[n2]
            L = np.sqrt((x2-x1)**2 + (y2-y1)**2)
            c = (x2-x1) / L
            s = (y2-y1) / L
            
            # 单元刚度矩阵
            k = self.E * areas[i] / L * np.array([
                [c*c, c*s, -c*c, -c*s],
                [c*s, s*s, -c*s, -s*s],
                [-c*c, -c*s, c*c, c*s],
                [-c*s, -s*s, c*s, s*s]
            ])
            
            # 组装到全局矩阵
            dofs = [2*n1, 2*n1+1, 2*n2, 2*n2+1]
            for ii, di in enumerate(dofs):
                for jj, dj in enumerate(dofs):
                    K[di, dj] += k[ii, jj]
        
        # 求解
        free_dofs = np.setdiff1d(np.arange(self.ndof), fixed_dofs)
        U = np.zeros(self.ndof)
        
        try:
            U[free_dofs] = np.linalg.solve(K[free_dofs][:, free_dofs], loads[free_dofs])
        except:
            U[free_dofs] = 0
        
        # 计算应力和重量
        stress = np.zeros(self.n_elements)
        weight = 0
        
        for i, (n1, n2) in enumerate(self.elements):
            x1, y1 = self.nodes[n1]
            x2, y2 = self.nodes[n2]
            L = np.sqrt((x2-x1)**2 + (y2-y1)**2)
            c = (x2-x1) / L
            s = (y2-y1) / L
            
            dofs = [2*n1, 2*n1+1, 2*n2, 2*n2+1]
            u = U[dofs]
            
            # 应变 = B * u
            B = np.array([-c, -s, c, s]) / L
            strain = B @ u
            stress[i] = self.E * strain
            
            weight += self.rho * areas[i] * L
        
        return U, stress, weight


class SurrogateModel:
    """代理模型（神经网络）"""
    
    def __init__(self, input_dim, output_dim, hidden_dims=[64, 64]):
        self.input_dim = input_dim
        self.output_dim = output_dim
        self.hidden_dims = hidden_dims
        
        # 初始化权重
        self.weights = []
        self.biases = []
        
        dims = [input_dim] + hidden_dims + [output_dim]
        for i in range(len(dims) - 1):
            self.weights.append(np.random.randn(dims[i], dims[i+1]) * 0.1)
            self.biases.append(np.zeros((1, dims[i+1])))
        
        # 存储训练数据用于归一化
        self.X_mean = None
        self.X_std = None
        self.y_mean = None
        self.y_std = None
    
    def relu(self, x):
        return np.maximum(0, x)
    
    def forward(self, X):
        """前向传播"""
        a = X
        for i, (W, b) in enumerate(zip(self.weights, self.biases)):
            z = a @ W + b
            if i < len(self.weights) - 1:
                a = self.relu(z)
            else:
                a = z
        return a
    
    def train(self, X, y, epochs=500, lr=0.001):
        """训练模型"""
        # 数据归一化
        self.X_mean = np.mean(X, axis=0)
        self.X_std = np.std(X, axis=0) + 1e-8
        self.y_mean = np.mean(y, axis=0)
        self.y_std = np.std(y, axis=0) + 1e-8
        
        X_norm = (X - self.X_mean) / self.X_std
        y_norm = (y - self.y_mean) / self.y_std
        
        losses = []
        
        for epoch in range(epochs):
            # 前向传播
            activations = [X_norm]
            a = X_norm
            for i, (W, b) in enumerate(zip(self.weights, self.biases)):
                z = a @ W + b
                if i < len(self.weights) - 1:
                    a = self.relu(z)
                else:
                    a = z
                activations.append(a)
            
            # 计算损失
            loss = np.mean((a - y_norm)**2)
            losses.append(loss)
            
            # 反向传播（简化版）
            delta = 2 * (a - y_norm) / len(X_norm)
            
            for i in range(len(self.weights) - 1, -1, -1):
                dW = activations[i].T @ delta
                db = np.sum(delta, axis=0, keepdims=True)
                
                self.weights[i] -= lr * dW
                self.biases[i] -= lr * db
                
                if i > 0:
                    delta = delta @ self.weights[i].T
                    delta *= (activations[i] > 0).astype(float)
            
            if (epoch + 1) % 100 == 0:
                print(f"  Epoch {epoch+1}/{epochs}, Loss: {loss:.6f}")
        
        return losses
    
    def predict(self, X):
        """预测"""
        X_norm = (X - self.X_mean) / self.X_std
        y_norm = self.forward(X_norm)
        return y_norm * self.y_std + self.y_mean


def create_10_bar_truss():
    """创建10杆桁架结构"""
    # 节点坐标
    nodes = np.array([
        [0, 0],      # 0
        [360, 0],    # 1
        [0, 360],    # 2
        [360, 360],  # 3
        [720, 360],  # 4
        [720, 0]     # 5
    ]) * 0.0254  # 转换为米
    
    # 单元连接
    elements = np.array([
        [0, 1], [1, 2], [2, 3], [3, 0],  # 底部和垂直
        [0, 2], [1, 3],                  # 对角线
        [1, 4], [3, 4], [4, 5], [3, 5]   # 右侧
    ])
    
    return nodes, elements


def generate_surrogate_data(fem, n_samples=500):
    """生成代理模型训练数据"""
    print("生成代理模型训练数据...")
    
    # 载荷条件
    loads = np.zeros(fem.ndof)
    loads[2*4+1] = -1e5  # 节点4y方向向下
    loads[2*5+1] = -1e5  # 节点5y方向向下
    
    # 固定节点0和5
    fixed_dofs = [0, 1, 2*5, 2*5+1]
    
    X = []
    y = []
    
    # 随机采样截面积
    np.random.seed(42)
    for i in range(n_samples):
        areas = np.random.uniform(1e-4, 5e-4, fem.n_elements)
        
        U, stress, weight = fem.analyze(areas, loads, fixed_dofs)
        
        max_disp = np.max(np.abs(U))
        max_stress = np.max(np.abs(stress))
        
        X.append(areas)
        y.append([weight, max_disp, max_stress])
        
        if (i + 1) % 100 == 0:
            print(f"  已生成 {i+1}/{n_samples} 个样本")
    
    return np.array(X), np.array(y)


def optimize_with_surrogate(fem, surrogate, max_iter=50):
    """使用代理模型进行优化"""
    print("\n使用代理模型进行优化...")
    
    # 初始化
    areas = np.ones(fem.n_elements) * 2e-4
    
    # 约束
    max_disp_allowed = 0.05
    max_stress_allowed = 2.5e8
    
    history = {'weight': [], 'max_disp': [], 'max_stress': []}
    
    lr = 1e-6
    
    for iter in range(max_iter):
        # 使用代理模型预测
        pred = surrogate.predict(areas.reshape(1, -1))[0]
        weight, max_disp, max_stress = pred
        
        # 数值梯度
        grad = np.zeros(fem.n_elements)
        eps = 1e-6
        
        for i in range(fem.n_elements):
            areas_plus = areas.copy()
            areas_plus[i] += eps
            pred_plus = surrogate.predict(areas_plus.reshape(1, -1))[0]
            
            # 目标：最小化重量，惩罚约束违反
            penalty = 0
            if pred_plus[1] > max_disp_allowed:
                penalty += 1e6 * (pred_plus[1] - max_disp_allowed)**2
            if pred_plus[2] > max_stress_allowed:
                penalty += 1e6 * (pred_plus[2] - max_stress_allowed)**2
            
            obj_plus = pred_plus[0] + penalty
            obj = weight + (1e6 * (max(max_disp - max_disp_allowed, 0)**2 + 
                                   max(max_stress - max_stress_allowed, 0)**2))
            
            grad[i] = (obj_plus - obj) / eps
        
        # 更新
        areas -= lr * grad
        areas = np.clip(areas, 1e-4, 5e-4)
        
        history['weight'].append(weight)
        history['max_disp'].append(max_disp)
        history['max_stress'].append(max_stress)
        
        if (iter + 1) % 10 == 0:
            print(f"  Iter {iter+1}: Weight={weight:.2f}kg, "
                  f"Disp={max_disp*1000:.2f}mm, Stress={max_stress/1e6:.2f}MPa")
    
    return areas, history


def optimize_traditional(fem, max_iter=50):
    """传统优化方法"""
    print("\n使用传统方法进行优化...")
    
    # 载荷和边界条件
    loads = np.zeros(fem.ndof)
    loads[2*4+1] = -1e5
    loads[2*5+1] = -1e5
    fixed_dofs = [0, 1, 2*5, 2*5+1]
    
    areas = np.ones(fem.n_elements) * 2e-4
    
    max_disp_allowed = 0.05
    max_stress_allowed = 2.5e8
    
    history = {'weight': [], 'max_disp': [], 'max_stress': []}
    
    lr = 1e-6
    
    for iter in range(max_iter):
        U, stress, weight = fem.analyze(areas, loads, fixed_dofs)
        
        max_disp = np.max(np.abs(U))
        max_stress_val = np.max(np.abs(stress))
        
        # 数值梯度
        grad = np.zeros(fem.n_elements)
        eps = 1e-6
        
        for i in range(fem.n_elements):
            areas_plus = areas.copy()
            areas_plus[i] += eps
            _, stress_plus, weight_plus = fem.analyze(areas_plus, loads, fixed_dofs)
            
            max_disp_plus = np.max(np.abs(U))  # 简化，实际应该重新计算
            max_stress_plus = np.max(np.abs(stress_plus))
            
            penalty = 0
            if max_disp_plus > max_disp_allowed:
                penalty += 1e6 * (max_disp_plus - max_disp_allowed)**2
            if max_stress_plus > max_stress_allowed:
                penalty += 1e6 * (max_stress_plus - max_stress_allowed)**2
            
            obj_plus = weight_plus + penalty
            obj = weight + (1e6 * (max(max_disp - max_disp_allowed, 0)**2 + 
                                   max(max_stress_val - max_stress_allowed, 0)**2))
            
            grad[i] = (obj_plus - obj) / eps
        
        areas -= lr * grad
        areas = np.clip(areas, 1e-4, 5e-4)
        
        history['weight'].append(weight)
        history['max_disp'].append(max_disp)
        history['max_stress'].append(max_stress_val)
        
        if (iter + 1) % 10 == 0:
            print(f"  Iter {iter+1}: Weight={weight:.2f}kg, "
                  f"Disp={max_disp*1000:.2f}mm, Stress={max_stress_val/1e6:.2f}MPa")
    
    return areas, history


def main():
    """主程序"""
    print("="*60)
    print("基于代理模型的结构优化加速")
    print("="*60)
    
    # 创建桁架结构
    nodes, elements = create_10_bar_truss()
    fem = TrussFEM(nodes, elements)
    
    print(f"桁架结构: {fem.n_nodes}个节点, {fem.n_elements}个单元")
    
    # 生成训练数据
    X_train, y_train = generate_surrogate_data(fem, n_samples=500)
    
    # 训练代理模型
    print("\n训练代理模型...")
    surrogate = SurrogateModel(input_dim=fem.n_elements, output_dim=3, 
                               hidden_dims=[64, 64, 32])
    losses = surrogate.train(X_train, y_train, epochs=500, lr=0.001)
    
    # 验证代理模型精度
    print("\n验证代理模型...")
    X_test, y_test = generate_surrogate_data(fem, n_samples=100)
    y_pred = surrogate.predict(X_test)
    
    mse = np.mean((y_test - y_pred)**2, axis=0)
    print(f"测试集MSE: Weight={mse[0]:.2e}, Disp={mse[1]:.2e}, Stress={mse[2]:.2e}")
    
    # 使用代理模型优化
    areas_surrogate, history_surrogate = optimize_with_surrogate(fem, surrogate, max_iter=50)
    
    # 传统优化
    areas_traditional, history_traditional = optimize_traditional(fem, max_iter=50)
    
    # 可视化
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 代理模型训练损失
    ax = axes[0, 0]
    ax.plot(losses, linewidth=2)
    ax.set_xlabel('Epoch')
    ax.set_ylabel('Loss')
    ax.set_title('代理模型训练损失')
    ax.grid(True, alpha=0.3)
    ax.set_yscale('log')
    
    # 重量对比
    ax = axes[0, 1]
    ax.plot(history_surrogate['weight'], 'b-', linewidth=2, label='代理模型优化')
    ax.plot(history_traditional['weight'], 'r--', linewidth=2, label='传统优化')
    ax.set_xlabel('Iteration')
    ax.set_ylabel('Weight (kg)')
    ax.set_title('结构重量优化历史')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 最大位移
    ax = axes[1, 0]
    ax.plot(history_surrogate['max_disp'], 'b-', linewidth=2, label='代理模型优化')
    ax.plot(history_traditional['max_disp'], 'r--', linewidth=2, label='传统优化')
    ax.axhline(y=0.05, color='g', linestyle=':', label='约束限值')
    ax.set_xlabel('Iteration')
    ax.set_ylabel('Max Displacement (m)')
    ax.set_title('最大位移优化历史')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # 最大应力
    ax = axes[1, 1]
    ax.plot(history_surrogate['max_stress'], 'b-', linewidth=2, label='代理模型优化')
    ax.plot(history_traditional['max_stress'], 'r--', linewidth=2, label='传统优化')
    ax.axhline(y=2.5e8, color='g', linestyle=':', label='约束限值')
    ax.set_xlabel('Iteration')
    ax.set_ylabel('Max Stress (Pa)')
    ax.set_title('最大应力优化历史')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('surrogate_optimization.png', dpi=150, bbox_inches='tight')
    print("\n结果图已保存为 surrogate_optimization.png")
    
    print("\n" + "="*60)
    print("程序执行完成！")
    print("="*60)


if __name__ == "__main__":
    main()

3.2.3 代码解析

桁架有限元分析：

class TrussFEM:
    def analyze(self, areas, loads, fixed_dofs):
        # 组装刚度矩阵
        K = np.zeros((self.ndof, self.ndof))
        
        for i, (n1, n2) in enumerate(self.elements):
            x1, y1 = self.nodes[n1]
            x2, y2 = self.nodes[n2]
            L = np.sqrt((x2-x1)**2 + (y2-y1)**2)
            c = (x2-x1) / L
            s = (y2-y1) / L
            
            k = self.E * areas[i] / L * np.array([...])
            
            dofs = [2*n1, 2*n1+1, 2*n2, 2*n2+1]
            for ii, di in enumerate(dofs):
                for jj, dj in enumerate(dofs):
                    K[di, dj] += k[ii, jj]

这段代码实现了桁架结构的有限元分析。对于每个单元，计算其长度、方向余弦，然后组装单元刚度矩阵到全局刚度矩阵。

代理模型训练：

def generate_surrogate_data(fem, n_samples=500):
    X = []
    y = []
    
    for i in range(n_samples):
        areas = np.random.uniform(1e-4, 5e-4, fem.n_elements)
        U, stress, weight = fem.analyze(areas, loads, fixed_dofs)
        
        max_disp = np.max(np.abs(U))
        max_stress = np.max(np.abs(stress))
        
        X.append(areas)
        y.append([weight, max_disp, max_stress])

通过随机采样不同的截面积组合，使用有限元分析计算对应的结构响应，生成代理模型的训练数据。

基于代理模型的优化：

def optimize_with_surrogate(fem, surrogate, max_iter=50):
    for iter in range(max_iter):
        # 使用代理模型预测
        pred = surrogate.predict(areas.reshape(1, -1))[0]
        weight, max_disp, max_stress = pred
        
        # 数值梯度
        grad = np.zeros(fem.n_elements)
        for i in range(fem.n_elements):
            areas_plus = areas.copy()
            areas_plus[i] += eps
            pred_plus = surrogate.predict(areas_plus.reshape(1, -1))[0]
            grad[i] = (obj_plus - obj) / eps
        
        areas -= lr * grad

使用代理模型替代有限元分析，通过数值梯度法进行优化。由于代理模型的推理速度远快于有限元分析，整个优化过程可以大幅加速。

3.2.4 运行结果预期

运行代码后，将生成以下结果：

1. 代理模型训练损失曲线：显示训练过程中损失值的下降趋势
2. 结构重量优化历史：对比代理模型优化和传统优化的收敛过程
3. 最大位移优化历史：显示约束条件的满足情况
4. 最大应力优化历史：显示应力约束的满足情况

预期代理模型优化方法能够达到与传统方法相近的优化结果，但计算时间显著减少。

3.3 案例三：物理信息神经网络（PINN）在结构分析中的应用

本案例展示如何使用PINN方法求解简单的结构力学问题。

3.3.1 问题描述

考虑一维弹性杆的轴向变形问题。杆的长度为 $L$，弹性模量为 $E$，截面积为 $A$，承受轴向分布载荷 $q(x)$。控制方程为：

$$EA\frac{d^{2}u}{d{x}^2}+q(x)=0$$

边界条件：$u(0)=0$（固定端），$EA\frac{du}{dx}{|}_{x=L}=F$（自由端受力）

3.3.2 完整代码实现

import numpy as np
import matplotlib.pyplot as plt

plt.switch_backend('Agg')


class PINN_1D_Bar:
    """
    物理信息神经网络求解一维杆件问题
    """
    
    def __init__(self, L=1.0, E=1.0, A=1.0, hidden_layers=[20, 20, 20]):
        """
        初始化PINN
        
        参数:
            L: 杆长
            E: 弹性模量
            A: 截面积
            hidden_layers: 隐藏层神经元数量
        """
        self.L = L
        self.E = E
        self.A = A
        self.EA = E * A
        
        # 网络结构: 1输入 -> 隐藏层 -> 1输出
        layers = [1] + hidden_layers + [1]
        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 = []
    
    def tanh(self, x):
        """双曲正切激活函数"""
        return np.tanh(x)
    
    def tanh_derivative(self, x):
        """tanh导数"""
        return 1 - np.tanh(x)**2
    
    def forward(self, x):
        """
        前向传播
        
        参数:
            x: 位置坐标 (n_points, 1)
        
        返回:
            u: 位移预测
            cache: 缓存用于反向传播
        """
        cache = []
        a = x
        
        for i, (W, b) in enumerate(zip(self.weights, self.biases)):
            z = a @ W + b
            cache.append({'z': z, 'a': a})
            
            if i < len(self.weights) - 1:
                a = self.tanh(z)
            else:
                a = z  # 输出层无激活函数
        
        return a, cache
    
    def compute_derivatives(self, x, cache):
        """
        计算位移对位置的一阶和二阶导数
        
        使用自动微分思想
        """
        n_points = x.shape[0]
        
        # 计算 du/dx
        # 从输出层开始反向传播梯度
        du_dz = [np.ones((n_points, 1))]  # 输出层对自身的梯度
        
        for i in range(len(self.weights) - 1, 0, -1):
            du_da = du_dz[-1] @ self.weights[i].T
            dz_da_prev = self.tanh_derivative(cache[i-1]['z'])
            du_dz.append(du_da * dz_da_prev)
        
        du_dz.reverse()
        
        # du/dx = du/dz[0] @ W[0].T
        du_dx = du_dz[0] @ self.weights[0].T
        
        # 计算 d²u/dx²（数值微分近似）
        eps = 1e-5
        x_plus = x + eps
        x_minus = x - eps
        
        u_plus, _ = self.forward(x_plus)
        u_minus, _ = self.forward(x_minus)
        d2u_dx2 = (u_plus - 2 * cache[-1]['a'] + u_minus) / (eps**2)
        
        return du_dx, d2u_dx2
    
    def compute_loss(self, x_collocation, x_bc, u_bc, x_force, force_value, q_func):
        """
        计算PINN损失函数
        
        包含三部分:
        1. 控制方程残差
        2. 位移边界条件
        3. 力边界条件
        """
        # 前向传播
        u_pred, cache = self.forward(x_collocation)
        
        # 计算导数
        du_dx, d2u_dx2 = self.compute_derivatives(x_collocation, cache)
        
        # 1. 控制方程残差: EA * d²u/dx² + q(x) = 0
        q_vals = q_func(x_collocation)
        residual_pde = self.EA * d2u_dx2 + q_vals
        loss_pde = np.mean(residual_pde**2)
        
        # 2. 位移边界条件: u(0) = 0
        u_bc_pred, _ = self.forward(x_bc)
        loss_bc = np.mean((u_bc_pred - u_bc)**2)
        
        # 3. 力边界条件: EA * du/dx = F at x=L
        du_dx_force, _ = self.compute_derivatives(x_force, cache)
        residual_force = self.EA * du_dx_force - force_value
        loss_force = np.mean(residual_force**2)
        
        # 总损失
        loss = loss_pde + 10.0 * loss_bc + 10.0 * loss_force
        
        return loss, {'pde': loss_pde, 'bc': loss_bc, 'force': loss_force}
    
    def train(self, n_collocation=100, epochs=5000, lr=0.001):
        """
        训练PINN
        
        参数:
            n_collocation: 配点数量
            epochs: 训练轮数
            lr: 学习率
        """
        # 生成配点
        x_collocation = np.random.uniform(0, self.L, (n_collocation, 1))
        
        # 边界条件点
        x_bc = np.array([[0.0]])  # 固定端
        u_bc = np.array([[0.0]])
        
        x_force = np.array([[self.L]])  # 自由端
        force_value = np.array([[1.0]])  # 单位力
        
        # 分布载荷函数
        q_func = lambda x: np.zeros_like(x)
        
        print(f"开始训练PINN，共{epochs}轮...")
        
        for epoch in range(epochs):
            # 计算损失
            loss, loss_dict = self.compute_loss(
                x_collocation, x_bc, u_bc, x_force, force_value, q_func
            )
            self.loss_history.append(loss)
            
            # 数值梯度下降（简化版）
            eps = 1e-5
            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_collocation, x_bc, u_bc, x_force, force_value, q_func
                        )
                        grad = (loss_plus - loss) / eps
                        self.weights[i][j, k] -= eps
                        self.weights[i][j, k] -= lr * grad
            
            if (epoch + 1) % 500 == 0:
                print(f"  Epoch {epoch+1}: Loss={loss:.6f}, "
                      f"PDE={loss_dict['pde']:.6f}, BC={loss_dict['bc']:.6f}")
        
        print("训练完成！")
    
    def predict(self, x):
        """预测位移"""
        u, _ = self.forward(x)
        return u


def analytical_solution(x, L, E, A, F):
    """解析解"""
    return F * x / (E * A)


def main():
    """主程序"""
    print("="*60)
    print("物理信息神经网络（PINN）求解一维杆件问题")
    print("="*60)
    
    # 问题参数
    L = 1.0  # 杆长
    E = 1.0  # 弹性模量
    A = 1.0  # 截面积
    F = 1.0  # 端部力
    
    # 创建PINN
    pinn = PINN_1D_Bar(L=L, E=E, A=A, hidden_layers=[20, 20, 20])
    
    # 训练
    pinn.train(n_collocation=100, epochs=2000, lr=0.01)
    
    # 预测
    x_test = np.linspace(0, L, 100).reshape(-1, 1)
    u_pinn = pinn.predict(x_test)
    u_exact = analytical_solution(x_test, L, E, A, F)
    
    # 计算误差
    mse = np.mean((u_pinn - u_exact)**2)
    print(f"\n预测MSE: {mse:.6e}")
    
    # 可视化
    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    
    # 损失历史
    ax = axes[0]
    ax.plot(pinn.loss_history, linewidth=2)
    ax.set_xlabel('Epoch')
    ax.set_ylabel('Loss')
    ax.set_title('PINN训练损失')
    ax.grid(True, alpha=0.3)
    ax.set_yscale('log')
    
    # 位移对比
    ax = axes[1]
    ax.plot(x_test, u_exact, 'b-', linewidth=2, label='解析解')
    ax.plot(x_test, u_pinn, 'r--', linewidth=2, label='PINN预测')
    ax.set_xlabel('位置 x')
    ax.set_ylabel('位移 u')
    ax.set_title('位移分布对比')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig('pinn_solution.png', dpi=150, bbox_inches='tight')
    print("结果图已保存为 pinn_solution.png")
    
    print("\n" + "="*60)
    print("程序执行完成！")
    print("="*60)


if __name__ == "__main__":
    main()

3.3.3 代码解析

PINN损失函数：

def compute_loss(self, x_collocation, x_bc, u_bc, x_force, force_value, q_func):
    # 1. 控制方程残差
    residual_pde = self.EA * d2u_dx2 + q_vals
    loss_pde = np.mean(residual_pde**2)
    
    # 2. 位移边界条件
    loss_bc = np.mean((u_bc_pred - u_bc)**2)
    
    # 3. 力边界条件
    loss_force = np.mean(residual_force**2)
    
    # 总损失
    loss = loss_pde + 10.0 * loss_bc + 10.0 * loss_force

PINN的核心思想是将物理约束嵌入损失函数。这里包括三部分损失：控制方程残差（PDE）、位移边界条件（BC）和力边界条件。通过调整权重系数，可以控制各项损失的相对重要性。

3.3.4 运行结果预期

运行代码后，将生成PINN训练损失曲线和位移对比图，验证PINN在求解结构力学问题中的有效性。