主题004：优化问题的数学表述

摘要

本教程系统介绍结构优化问题的数学表述方法，这是将工程优化问题转化为可计算数学模型的关键步骤。内容涵盖优化问题的标准形式、设计变量与约束条件的定义、目标函数的构建方法。详细讲解等式约束与不等式约束的处理技巧，以及约束函数的归一化方法。介绍灵敏度分析的基本概念及其在优化中的重要性。通过Python实现，演示如何将实际工程问题（如桁架重量最小化）转化为标准优化问题，并使用SciPy进行求解。本教程为后续各类结构优化算法的学习奠定数学基础。

关键词

优化问题；数学表述；设计变量；约束条件；目标函数；灵敏度分析；归一化；Python优化

---

1. 引言

1.1 为什么需要数学表述

工程优化问题往往以自然语言或物理描述的形式出现，例如：

• “设计一个重量最轻的桁架，使其能够承受给定载荷”
• “优化梁的截面尺寸，在满足强度要求的前提下降低成本”
• “确定最佳的结构形状，使应力分布均匀”

要将这些问题转化为计算机可以求解的形式，必须建立严格的数学模型。数学表述是连接工程需求与数值算法的桥梁。

1.2 优化问题的基本要素

任何优化问题都包含三个核心要素：

1. 设计变量（Design Variables）：可以调整的参数，如尺寸、形状、材料分布等
2. 目标函数（Objective Function）：需要最小化或最大化的性能指标
3. 约束条件（Constraints）：必须满足的限制条件，如强度、刚度、稳定性等

1.3 本教程的学习目标

通过本教程的学习，读者将能够：

1. 掌握优化问题的标准数学形式：理解min-max问题的表述规范

2. 正确选择设计变量：根据问题特点确定合适的优化参数

3. 构建合理的目标函数：将工程需求转化为数学表达式

4. 处理各类约束条件：掌握等式约束、不等式约束的数学表述

5. 实现优化问题的Python建模：使用SciPy等工具求解标准优化问题

---

2. 优化问题的标准形式

2.1 一般数学表述

结构优化问题的标准数学形式为：

$$\begin{array}{rlrl}& \underset{\mathbf{x}}{\mathrm{minimize}}& & f(\mathbf{x})\\ & \text{subject to}& & g_j(\mathbf{x})\le 0,j=1,2,\dots ,m\\ & & & h_k(\mathbf{x})=0,k=1,2,\dots ,l\\ & & & x_i^L\le x_i\le x_i^U,i=1,2,\dots ,n\end{array}$$

其中：

• $\mathbf{x}=[x_1,x_2,\dots ,x_n{]}^T$：设计变量向量
• $f(\mathbf{x})$：目标函数
• $g_j(\mathbf{x})$：不等式约束函数
• $h_k(\mathbf{x})$：等式约束函数
• $x_i^L,x_i^U$：设计变量的上下界

2.2 设计空间与可行域

设计空间（Design Space）：由所有设计变量构成的$n$维空间，每个点代表一个设计方案。

可行域（Feasible Region）：满足所有约束条件的设计点集合，记为：

F={x∈Rn∣gj(x)≤0,hk(x)=0,xiL≤xi≤xiU} \mathcal{F} = \{\mathbf{x} \in \mathbb{R}^n \mid g_j(\mathbf{x}) \leq 0, h_k(\mathbf{x}) = 0, x_i^L \leq x_i \leq x_i^U\} F={x∈Rn∣gj​(x)≤0,hk​(x)=0,xiL​≤xi​≤xiU​}

2.3 最优解的分类

局部最优解（Local Optimum）：在邻域内目标函数值最小的可行解

全局最优解（Global Optimum）：在整个可行域内目标函数值最小的解

凸优化问题：若目标函数和可行域都是凸的，则局部最优即为全局最优。

---

3. 设计变量的选择

3.1 设计变量的类型

3.1.1 连续变量

可以在给定区间内取任意实数值，如：

• 杆件截面面积 $A\in [A_{min},A_{max}]$
• 板件厚度 $t\in [t_{min},t_{max}]$
• 节点坐标 $(x,y,z)$

3.1.2 离散变量

只能取特定离散值，如：

• 标准型材规格（工字钢型号）
• 材料类型（钢、铝、复合材料）
• 加强筋数量（整数）

3.1.3 混合变量

同时包含连续和离散变量的优化问题，如拓扑优化中的0-1变量。

3.2 设计变量的尺度问题

无量纲化：当设计变量量纲和量级差异较大时，应进行归一化处理：

$${\tilde{x}}_i=\frac{x_i-x_i^L}{x_i^U-x_i^L}$$

归一化后的变量 ${\tilde{x}}_i\in [0,1]$，有利于优化算法的收敛。

3.3 设计变量的耦合

显式耦合：变量间存在显式关系，如 $x_1+x_2=C$

隐式耦合：通过约束条件或物理方程间接关联

---

4. 目标函数的构建

4.1 常见目标函数类型

4.1.1 重量最小化

$$f(\mathbf{x})=\sum _{e=1}^{n_e}{\rho }_e{V}_e(\mathbf{x})\to \min$$

其中 ${\rho }_e$ 是材料密度，$V_e$ 是单元体积。

4.1.2 刚度最大化（柔度最小化）

$$f(\mathbf{x})={\mathbf{F}}^T\mathbf{U}={\mathbf{U}}^T\mathbf{KU}\to \min$$

其中 $\mathbf{F}$ 是载荷向量，$\mathbf{U}$ 是位移向量，$\mathbf{K}$ 是刚度矩阵。

4.1.3 频率优化

最大化基频：

$$f(\mathbf{x})=-{\omega }_1(\mathbf{x})\to \min$$

或使频率避开特定区间。

4.1.4 多目标优化

当存在多个冲突目标时，采用加权求和或Pareto最优：

$$f(\mathbf{x})=w_1{f}_1(\mathbf{x})+w_2{f}_2(\mathbf{x})+\cdots +w_m{f}_m(\mathbf{x})$$

4.2 目标函数的尺度变换

归一化处理：将不同量纲的目标函数归一化到同一量级：

$$\tilde{f}(\mathbf{x})=\frac{f(\mathbf{x})-f_{ref}}{f_{ref}}$$

其中 $f_{ref}$ 是参考值（如初始设计的函数值）。

---

5. 约束条件的处理

5.1 约束类型

5.1.1 应力约束

$$g(\mathbf{x})=\frac{{\sigma }_{max}(\mathbf{x})}{[\sigma ]}-1\le 0$$

其中 $[\sigma ]$ 是许用应力。

5.1.2 位移约束

$$g(\mathbf{x})=\frac{|u_i(\mathbf{x})|}{[u]}-1\le 0$$

其中 $[u]$ 是许用位移。

5.1.3 稳定性约束

$$g(\mathbf{x})=1-\frac{P_{cr}(\mathbf{x})}{P_{applied}}\le 0$$

其中 $P_{cr}$ 是临界屈曲载荷。

5.1.4 几何约束

制造和装配要求：

$$g(\mathbf{x})=x_{min}-x_i\le 0\text{或}{x}_i-x_{max}\le 0$$

5.2 约束归一化

重要性：不同约束的量纲和量级差异很大，需要归一化处理。

归一化方法：

$$\tilde{g}(\mathbf{x})=\frac{g(\mathbf{x})}{g_{allow}}-1\le 0$$

归一化后，所有约束都在同一量级，便于优化算法处理。

5.3 约束的激活与违反

激活约束（Active Constraint）：$g_j(\mathbf{x})=0$，在最优解处通常存在激活约束

违反约束（Violated Constraint）：$g_j(\mathbf{x})>0$，不可行解

松弛约束（Inactive Constraint）：$g_j(\mathbf{x})<0$，不起限制作用

---

6. 灵敏度分析

6.1 灵敏度的定义

设计灵敏度：目标函数或约束函数对设计变量的导数：

$$\frac{\partial f}{\partial x_i},\frac{\partial g_j}{\partial x_i}$$

6.2 灵敏度的计算方法

6.2.1 有限差分法

$$\frac{\partial f}{\partial x_i}\approx \frac{f(x_i+\Delta x_i)-f(x_i)}{\Delta x_i}$$

优点：实现简单；缺点：计算量大，精度受步长影响。

6.2.2 解析法

利用链式法则直接推导解析表达式，精度高但推导复杂。

6.2.3 伴随法（Adjoint Method）

对于大量设计变量的情况，伴随法效率最高：

$$\frac{\partial f}{\partial \mathbf{x}}=\frac{\partial f}{\partial \mathbf{u}}\frac{\partial \mathbf{u}}{\partial \mathbf{x}}$$

6.3 灵敏度在优化中的作用

1. 指导搜索方向：梯度信息用于确定优化方向
2. 判断约束重要性：灵敏度大的约束对设计影响更大
3. 近似分析：用于构建近似模型，减少重分析次数

---

7. Python实现：桁架重量优化

7.1 问题描述

考虑一个10杆平面桁架结构，优化目标是最小化结构重量，约束条件包括：

• 各杆件应力不超过许用应力
• 节点位移不超过许用位移
• 杆件截面积有上下限

7.2 完整Python代码

# -*- coding: utf-8 -*-
"""
主题004：优化问题的数学表述
案例：10杆桁架重量优化
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.optimize import minimize
import warnings
warnings.filterwarnings('ignore')
import os

# 创建输出目录
output_dir = r'd:\文档\500仿真领域\工程仿真\结构优化设计\主题004'
os.makedirs(output_dir, exist_ok=True)

plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 60)
print("优化问题的数学表述 - 10杆桁架重量优化")
print("=" * 60)

# ==================== 1. 问题参数定义 ====================
print("\n【步骤1】定义优化问题参数...")

# 材料参数
E = 1e7       # 弹性模量，psi
rho = 0.1     # 密度，lb/in³

# 许用应力（拉压相同）
sigma_allow = 25000  # psi

# 许用位移
u_allow = 2.0  # in

# 几何参数
L = 360.0  # in，水平杆件长度
H = 360.0  # in，竖直方向高度

print(f"  弹性模量 E = {E/1e6:.1f} × 10⁶ psi")
print(f"  材料密度 ρ = {rho} lb/in³")
print(f"  许用应力 [σ] = {sigma_allow/1000:.1f} ksi")
print(f"  许用位移 [u] = {u_allow} in")
print(f"  杆件长度 L = {L} in")

# ==================== 2. 定义有限元分析函数 ====================
print("\n【步骤2】定义有限元分析函数...")

def truss_10bar_fea(A):
    """
    10杆桁架有限元分析
    返回：节点位移、杆件应力
    """
    # 节点坐标 (6个节点，2维)
    nodes = np.array([
        [2*L, H],    # 节点1
        [L, H],      # 节点2
        [0, H],      # 节点3
        [2*L, 0],    # 节点4
        [L, 0],      # 节点5
        [0, 0]       # 节点6
    ])
    
    # 杆件连接关系 (10根杆)
    elements = np.array([
        [0, 1],   # 杆1
        [1, 2],   # 杆2
        [2, 5],   # 杆3
        [5, 4],   # 杆4
        [4, 3],   # 杆5
        [3, 0],   # 杆6
        [1, 4],   # 杆7
        [2, 4],   # 杆8
        [1, 5],   # 杆9
        [0, 4]    # 杆10
    ])
    
    n_nodes = len(nodes)
    n_elements = len(elements)
    n_dofs = 2 * n_nodes
    
    # 组装整体刚度矩阵
    K = np.zeros((n_dofs, n_dofs))
    
    for e in range(n_elements):
        node_i, node_j = elements[e]
        xi, yi = nodes[node_i]
        xj, yj = nodes[node_j]
        
        L_e = np.sqrt((xj - xi)**2 + (yj - yi)**2)
        c = (xj - xi) / L_e
        s = (yj - yi) / L_e
        
        # 单元刚度矩阵（局部坐标）
        k_local = (E * A[e] / L_e) * np.array([[1, -1], [-1, 1]])
        
        # 转换矩阵
        T = np.array([[c, s, 0, 0], [0, 0, c, s]])
        
        # 单元刚度矩阵（全局坐标）
        k_global = T.T @ k_local @ T
        
        # 自由度映射
        dof_i = [2*node_i, 2*node_i+1]
        dof_j = [2*node_j, 2*node_j+1]
        dof_map = dof_i + dof_j
        
        # 组装
        for i in range(4):
            for j in range(4):
                K[dof_map[i], dof_map[j]] += k_global[i, j]
    
    # 载荷向量
    F = np.zeros(n_dofs)
    F[2*1+1] = -1e5  # 节点2的y方向载荷
    F[2*3+1] = -1e5  # 节点4的y方向载荷
    
    # 边界条件（固定节点5和6）
    fixed_dofs = [2*4, 2*4+1, 2*5, 2*5+1]  # 节点5和6的所有自由度
    free_dofs = [i for i in range(n_dofs) if i not in fixed_dofs]
    
    # 求解位移
    K_ff = K[np.ix_(free_dofs, free_dofs)]
    F_f = F[free_dofs]
    
    try:
        U_f = np.linalg.solve(K_ff, F_f)
    except:
        return None, None
    
    U = np.zeros(n_dofs)
    U[free_dofs] = U_f
    
    # 计算杆件应力
    stresses = np.zeros(n_elements)
    for e in range(n_elements):
        node_i, node_j = elements[e]
        xi, yi = nodes[node_i]
        xj, yj = nodes[node_j]
        
        L_e = np.sqrt((xj - xi)**2 + (yj - yi)**2)
        c = (xj - xi) / L_e
        s = (yj - yi) / L_e
        
        # 提取单元位移
        u_e = np.array([U[2*node_i], U[2*node_i+1], U[2*node_j], U[2*node_j+1]])
        
        # 应变 = (u_j - u_i) / L
        epsilon = np.array([-c, -s, c, s]) @ u_e / L_e
        
        # 应力
        stresses[e] = E * epsilon
    
    return U, stresses

# 测试有限元分析
A_test = np.ones(10) * 10.0  # 初始截面积
U_test, stress_test = truss_10bar_fea(A_test)

if U_test is not None:
    print(f"  ✓ 有限元分析函数测试成功")
    print(f"  初始设计最大位移：{np.max(np.abs(U_test)):.4f} in")
    print(f"  初始设计最大应力：{np.max(np.abs(stress_test))/1000:.2f} ksi")
else:
    print(f"  ✗ 有限元分析失败")

# ==================== 3. 定义优化问题 ====================
print("\n【步骤3】定义优化问题...")

# 设计变量：10根杆的截面积
# 约束：应力约束（20个，每根杆拉压各一个）
#       位移约束（4个，节点1和3的x,y方向）

# 设计变量边界
A_min = 0.1   # in²
A_max = 50.0  # in²
bounds = [(A_min, A_max) for _ in range(10)]

print(f"  设计变量数：10（杆件截面积）")
print(f"  截面积范围：[{A_min}, {A_max}] in²")
print(f"  不等式约束数：24（20个应力约束 + 4个位移约束）")

def objective(A):
    """目标函数：结构总重量"""
    # 杆件长度
    lengths = np.array([
        L, L, L, L, L, L,           # 杆1-6
        np.sqrt(2)*L, np.sqrt(2)*L, # 杆7-8
        np.sqrt(2)*L, np.sqrt(2)*L  # 杆9-10
    ])
    
    # 总重量 = 密度 × 截面积 × 长度
    weight = rho * np.sum(A * lengths)
    return weight

def stress_constraints(A):
    """应力约束：g_i = |σ_i|/[σ] - 1 ≤ 0"""
    U, stresses = truss_10bar_fea(A)
    
    if U is None:
        return np.ones(20) * 1e10  # 返回大值表示不可行
    
    # 应力约束（拉压各一个）
    g_stress = np.abs(stresses) / sigma_allow - 1.0
    
    return g_stress

def displacement_constraints(A):
    """位移约束：g_i = |u_i|/[u] - 1 ≤ 0"""
    U, stresses = truss_10bar_fea(A)
    
    if U is None:
        return np.ones(4) * 1e10
    
    # 提取节点1和3的位移
    u1x, u1y = U[0], U[1]   # 节点1
    u3x, u3y = U[4], U[5]   # 节点3
    
    # 位移约束
    g_disp = np.array([
        np.abs(u1x) / u_allow - 1.0,
        np.abs(u1y) / u_allow - 1.0,
        np.abs(u3x) / u_allow - 1.0,
        np.abs(u3y) / u_allow - 1.0
    ])
    
    return g_disp

def all_constraints(A):
    """所有约束函数"""
    g_stress = stress_constraints(A)
    g_disp = displacement_constraints(A)
    return np.concatenate([g_stress, g_disp])

# ==================== 4. 求解优化问题 ====================
print("\n【步骤4】使用SLSQP算法求解优化问题...")

# 初始设计
A0 = np.ones(10) * 10.0

# 约束条件（形式：g(x) ≤ 0）
constraints = {'type': 'ineq', 'fun': lambda A: -all_constraints(A)}

# 优化选项
options = {
    'ftol': 1e-6,
    'disp': True,
    'maxiter': 200
}

print(f"  初始重量：{objective(A0):.2f} lb")
print(f"  开始优化...")

# 求解优化问题
result = minimize(
    objective,
    A0,
    method='SLSQP',
    bounds=bounds,
    constraints=constraints,
    options=options
)

# ==================== 5. 结果分析 ====================
print("\n【步骤5】分析优化结果...")

if result.success:
    A_opt = result.x
    weight_opt = result.fun
    
    # 分析最优设计
    U_opt, stress_opt = truss_10bar_fea(A_opt)
    g_opt = all_constraints(A_opt)
    
    print(f"\n  ✓ 优化成功！")
    print(f"\n  【最优设计】")
    print(f"  最优重量：{weight_opt:.2f} lb")
    print(f"  重量减轻：{(1 - weight_opt/objective(A0))*100:.1f}%")
    
    print(f"\n  【最优截面积】")
    for i in range(10):
        print(f"  杆{i+1:2d}: A = {A_opt[i]:8.4f} in²")
    
    print(f"\n  【应力分析】")
    for i in range(10):
        print(f"  杆{i+1:2d}: σ = {stress_opt[i]/1000:8.2f} ksi " + 
              f"({'激活' if np.abs(g_opt[i]) < 0.01 else '松弛'})")
    
    print(f"\n  【位移分析】")
    print(f"  节点1 x方向: u = {U_opt[0]:8.4f} in")
    print(f"  节点1 y方向: u = {U_opt[1]:8.4f} in")
    print(f"  节点3 x方向: u = {U_opt[4]:8.4f} in")
    print(f"  节点3 y方向: u = {U_opt[5]:8.4f} in")
    
    print(f"\n  【约束违反情况】")
    max_violation = np.max(g_opt)
    if max_violation <= 0:
        print(f"  所有约束满足，最大约束值：{max_violation:.6f}")
    else:
        print(f"  警告：存在约束违反，最大违反：{max_violation:.6f}")
    
else:
    print(f"  ✗ 优化失败：{result.message}")

# ==================== 6. 可视化 ====================
print("\n【步骤6】生成可视化结果...")

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 子图1：截面积对比
ax1 = axes[0, 0]
x_pos = np.arange(1, 11)
width = 0.35
ax1.bar(x_pos - width/2, A0, width, label='Initial', alpha=0.7)
if result.success:
    ax1.bar(x_pos + width/2, A_opt, width, label='Optimal', alpha=0.7)
ax1.set_xlabel('Bar Number', fontsize=11)
ax1.set_ylabel('Cross-sectional Area (in²)', fontsize=11)
ax1.set_title('Cross-sectional Area Comparison', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 子图2：应力分布
ax2 = axes[0, 1]
if result.success:
    colors = ['red' if np.abs(g_opt[i]) < 0.01 else 'blue' for i in range(10)]
    ax2.bar(x_pos, np.abs(stress_opt)/1000, color=colors, alpha=0.7)
    ax2.axhline(y=sigma_allow/1000, color='r', linestyle='--', label='Allowable')
ax2.set_xlabel('Bar Number', fontsize=11)
ax2.set_ylabel('Stress (ksi)', fontsize=11)
ax2.set_title('Stress Distribution (Red = Active Constraint)', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)

# 子图3：优化迭代历史
ax3 = axes[1, 0]
if hasattr(result, 'nfev'):
    ax3.text(0.1, 0.8, f'Function evaluations: {result.nfev}', fontsize=12, transform=ax3.transAxes)
    ax3.text(0.1, 0.7, f'Iterations: {result.nit}', fontsize=12, transform=ax3.transAxes)
    ax3.text(0.1, 0.6, f'Initial weight: {objective(A0):.2f} lb', fontsize=12, transform=ax3.transAxes)
    if result.success:
        ax3.text(0.1, 0.5, f'Optimal weight: {result.fun:.2f} lb', fontsize=12, transform=ax3.transAxes)
        ax3.text(0.1, 0.4, f'Weight reduction: {(1 - result.fun/objective(A0))*100:.1f}%', fontsize=12, transform=ax3.transAxes)
ax3.set_xlim(0, 1)
ax3.set_ylim(0, 1)
ax3.axis('off')
ax3.set_title('Optimization Summary', fontsize=12, fontweight='bold')

# 子图4：约束违反情况
ax4 = axes[1, 1]
if result.success:
    # 应力约束
    g_stress = g_opt[:10]
    ax4.barh(x_pos, g_stress, alpha=0.7, label='Stress')
    ax4.axvline(x=0, color='r', linestyle='--')
    ax4.set_xlabel('Constraint Value', fontsize=11)
    ax4.set_ylabel('Bar Number', fontsize=11)
    ax4.set_title('Stress Constraint Violation', fontsize=12, fontweight='bold')
    ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/truss_optimization_results.png', dpi=150, bbox_inches='tight')
plt.close()
print(f"  ✓ 结果图已保存至: {output_dir}/truss_optimization_results.png")

# ==================== 7. 灵敏度分析 ====================
print("\n【步骤7】灵敏度分析...")

def compute_sensitivities(A, h=1e-4):
    """使用有限差分计算灵敏度"""
    n = len(A)
    
    # 目标函数灵敏度
    df_dA = np.zeros(n)
    f0 = objective(A)
    for i in range(n):
        A_perturbed = A.copy()
        A_perturbed[i] += h
        df_dA[i] = (objective(A_perturbed) - f0) / h
    
    return df_dA

if result.success:
    sensitivities = compute_sensitivities(A_opt)
    
    print(f"\n  目标函数对设计变量的灵敏度：")
    print(f"  {'杆件':<8} {'灵敏度':<15}")
    print(f"  {'-'*25}")
    for i in range(10):
        print(f"  {i+1:<8} {sensitivities[i]:<15.4f}")
    
    # 找出最敏感的设计变量
    max_sens_idx = np.argmax(np.abs(sensitivities))
    print(f"\n  最敏感的设计变量：杆{max_sens_idx+1}")
    print(f"  说明：改变该杆的截面积对总重量影响最大")

# ==================== 8. 结果总结 ====================
print("\n" + "=" * 60)
print("10杆桁架重量优化结果总结")
print("=" * 60)

print(f"\n【优化问题表述】")
print(f"  最小化：f(A) = Σ(ρ × A_i × L_i)")
print(f"  约束条件：")
print(f"    - 应力约束：|σ_i| ≤ {sigma_allow/1000:.1f} ksi")
print(f"    - 位移约束：|u_j| ≤ {u_allow} in")
print(f"    - 几何约束：{A_min} ≤ A_i ≤ {A_max} in²")

print(f"\n【求解方法】")
print(f"  算法：SLSQP (Sequential Least Squares Programming)")
print(f"  收敛准则：ftol = 1e-6")

if result.success:
    print(f"\n【优化结果】")
    print(f"  初始重量：{objective(A0):.2f} lb")
    print(f"  最优重量：{weight_opt:.2f} lb")
    print(f"  重量减轻：{(1 - weight_opt/objective(A0))*100:.1f}%")
    print(f"  迭代次数：{result.nit}")
    print(f"  函数评估次数：{result.nfev}")

print("=" * 60)

print("\n✓ 10杆桁架重量优化完成！")
print(f"✓ 所有结果已保存至: {output_dir}")

7.3 代码深度解析

7.3.1 有限元分析函数

truss_10bar_fea函数实现了10杆桁架的有限元分析，包括：

• 节点坐标和单元连接关系定义
• 整体刚度矩阵组装
• 边界条件处理
• 位移和应力计算

def truss_10bar_fea(A):
    """10杆桁架有限元分析"""
    # 节点坐标定义
    nodes = np.array([...])
    
    # 单元连接关系
    elements = np.array([...])
    
    # 组装刚度矩阵并求解
    ...
    return U, stresses

7.3.2 目标函数定义

目标函数计算结构总重量，是设计变量（截面积）的线性函数：

def objective(A):
    """目标函数：结构总重量"""
    lengths = np.array([...])  # 各杆长度
    weight = rho * np.sum(A * lengths)
    return weight

7.3.3 约束函数定义

约束函数采用归一化形式，便于优化算法统一处理：

def stress_constraints(A):
    """应力约束：g_i = |σ_i|/[σ] - 1 ≤ 0"""
    U, stresses = truss_10bar_fea(A)
    g_stress = np.abs(stresses) / sigma_allow - 1.0
    return g_stress

7.3.4 优化求解

使用SciPy的minimize函数，选择SLSQP算法：

result = minimize(
    objective,      # 目标函数
    A0,             # 初始值
    method='SLSQP', # 优化算法
    bounds=bounds,  # 变量边界
    constraints=constraints,  # 约束条件
    options=options # 优化选项
)

7.4 运行结果预期

运行代码后，将得到：

1. 优化前后对比：

   • 初始重量 vs 最优重量
   • 各杆截面积的变化
   • 重量减轻百分比

2. 约束分析：

   • 哪些约束被激活（达到边界）
   • 应力分布云图
   • 位移分布

3. 灵敏度信息：

   • 各设计变量对目标函数的影响程度
   • 指导进一步的设计改进

---