结构优化设计-主题005-优化算法基础

kkchenkx

90人浏览 · 2026-04-12 08:45:27

kkchenkx · 2026-04-12 08:45:27 发布

主题005：优化算法基础

一、引言：优化算法的核心地位

在结构优化设计中，优化算法是实现从"可行设计"到"最优设计"的关键桥梁。无论我们建立了多么精确的有限元模型，定义了多么合理的优化问题，如果没有高效的优化算法，所有的理论都将停留在纸面上。优化算法决定了我们能否在合理的时间内找到满足所有约束条件的最优解，以及这个最优解的质量如何。

优化算法的研究可以追溯到17世纪牛顿和莱布尼茨创立微积分时期，但真正的数值优化算法发展是在20世纪中叶随着计算机的出现而蓬勃兴起的。从早期的梯度下降法、牛顿法，到后来的共轭梯度法、拟牛顿法，再到近几十年兴起的遗传算法、粒子群优化、模拟退火等启发式算法，优化算法的大家族不断壮大，为解决各类复杂的工程优化问题提供了强大的工具箱。

在结构优化领域，优化算法的选择直接影响着设计效率和结果质量。尺寸优化问题通常具有明确的数学表达式和可计算的梯度信息，适合使用基于梯度的确定性算法；而拓扑优化问题往往涉及离散变量、非凸可行域，可能需要借助启发式算法或特殊的数学规划方法。理解不同优化算法的原理、特点和适用范围，是每一位结构优化工程师的必修课。

本教程将系统地介绍结构优化中常用的优化算法，从最基础的梯度下降法开始，逐步深入到牛顿法、共轭梯度法、遗传算法等经典算法。每个算法都将配有完整的Python实现和可视化演示，帮助读者直观地理解算法的工作原理和收敛行为。通过本教程的学习，读者将能够根据具体问题的特点选择合适的优化算法，并具备实现和调试优化算法的能力。
在这里插入图片描述

二、优化算法的分类体系

在深入具体算法之前，我们先建立一个清晰的算法分类框架。优化算法可以从多个维度进行分类：

2.1 按搜索策略分类

确定性算法（Deterministic Algorithms）：这类算法在相同的初始条件下总是产生相同的搜索路径和最终结果。它们通常基于数学分析，利用目标函数的梯度、Hessian矩阵等信息指导搜索方向。典型的确定性算法包括梯度下降法、牛顿法、共轭梯度法、拟牛顿法等。确定性算法的优点是收敛速度快、理论分析成熟，缺点是容易陷入局部最优，对初始点敏感。

随机性算法（Stochastic Algorithms）：这类算法在搜索过程中引入随机因素，即使从相同的初始点出发，也可能产生不同的搜索轨迹。随机性算法包括遗传算法、模拟退火、粒子群优化、蚁群算法等启发式方法。它们的优点是全局搜索能力强，能够跳出局部最优，适合处理非凸、多峰、离散优化问题；缺点是收敛速度慢，结果具有随机性，理论分析相对困难。

2.2 按梯度信息使用分类

基于梯度的算法（Gradient-based Methods）：这类算法需要计算目标函数和约束条件对设计变量的导数（梯度）。梯度的计算可以通过解析法（手工推导公式）、数值差分法或自动微分技术实现。基于梯度的算法收敛速度快，是光滑、连续优化问题的首选。在结构优化中，灵敏度分析（Sensitivity Analysis）就是计算结构响应（应力、位移等）对设计变量导数的过程。

无梯度算法（Gradient-free Methods）：这类算法仅利用函数值信息，不需要计算导数。它们适用于梯度难以计算或不存在的情况，如离散优化、黑箱优化、实验设计优化等。无梯度算法包括单纯形法、遗传算法、粒子群优化、Nelder-Mead方法等。在结构优化的早期探索阶段或概念设计阶段，无梯度算法往往更有优势。

2.3 按问题类型分类

无约束优化算法：处理没有任何约束条件的优化问题。虽然实际工程问题很少是完全无约束的，但许多约束优化算法通过惩罚函数法、拉格朗日乘子法等技术将约束问题转化为一系列无约束问题求解。

约束优化算法：直接处理带有等式约束和不等式约束的优化问题。序列二次规划（SQP）、内点法、可行方向法等都是经典的约束优化算法。

单目标优化算法：寻找使单一目标函数最优的设计方案。大多数经典优化算法都属于这一类。

多目标优化算法：同时优化多个相互冲突的目标函数，产生一组帕累托最优解供决策者选择。NSGA-II、MOEA/D等是多目标优化领域的代表性算法。

2.4 结构优化中的算法选择策略

在实际工程应用中，优化算法的选择需要综合考虑以下因素：

问题规模：设计变量的数量直接影响算法的选择。对于大规模问题（数千甚至数万个设计变量），基于梯度的算法通常是唯一可行的选择；对于小规模问题，启发式算法可能更合适。
函数光滑性：如果目标函数和约束条件都是光滑可微的，基于梯度的算法效率最高；如果存在不连续点、尖点或离散变量，则需要使用无梯度算法。
计算成本：结构分析（有限元计算）通常是优化过程中最耗时的环节。如果单次分析成本很高，应该优先选择收敛速度快的算法，减少分析次数。
全局最优需求：如果问题存在多个局部最优解，且需要找到全局最优，应该考虑使用全局优化算法或多起点策略。
约束复杂度：约束条件的数量、类型（线性/非线性、等式/不等式）也会影响算法选择。一些算法对约束处理有特殊优势，如SQP算法对非线性约束处理能力强。

三、梯度下降法：最基础的优化算法

3.1 算法原理与数学推导

梯度下降法（Gradient Descent）是最简单、最直观的优化算法，也是理解其他高级算法的基础。其核心思想是：函数在某点的梯度方向是函数值增长最快的方向，因此负梯度方向就是函数值下降最快的方向。

给定一个可微的目标函数 $f(x)f(\mathbf{x})$ ，其中 $x=[x1,x2,…,xn]T\mathbf{x} = [x_1, x_2, \ldots, x_n]^T$ 是设计变量向量。在点 $xk\mathbf{x}_k$ 处，函数的梯度为：

$∇f(xk)=[∂f∂x1,∂f∂x2,…,∂f∂xn]T\nabla f(\mathbf{x}_k) = \left[\frac{\partial f}{\partial x_1}, \frac{\partial f}{\partial x_2}, \ldots, \frac{\partial f}{\partial x_n}\right]^T$

梯度下降法的迭代公式为：

$xk+1=xk−αk∇f(xk)\mathbf{x}_{k+1} = \mathbf{x}_k - \alpha_k \nabla f(\mathbf{x}_k)$

其中， $αk>0\alpha_k > 0$ 是步长（学习率），控制着沿负梯度方向前进的距离。

3.2 步长选择策略

步长的选择对梯度下降法的性能至关重要：

固定步长：使用常数 $α\alpha$ ，实现简单但难以适应函数在不同区域的曲率变化。步长太大可能导致震荡或发散，步长太小则收敛缓慢。

精确线搜索：在每次迭代中，求解一维优化问题找到最优步长：
$αk=arg⁡min⁡α>0f(xk−α∇f(xk))\alpha_k = \arg\min_{\alpha > 0} f(\mathbf{x}_k - \alpha \nabla f(\mathbf{x}_k))$

精确线搜索保证每次迭代都沿搜索方向达到最小值，但计算成本较高。

回溯线搜索（Backtracking Line Search）：一种实用的近似线搜索策略。从一个较大的初始步长开始，以固定比例衰减（如 $α←βα\alpha \leftarrow \beta \alpha$ ，其中 $β∈(0,1)\beta \in (0,1)$ ），直到满足Armijo条件：

$f(xk−α∇f(xk))≤f(xk)−cα∥∇f(xk)∥2f(\mathbf{x}_k - \alpha \nabla f(\mathbf{x}_k)) \leq f(\mathbf{x}_k) - c \alpha \|\nabla f(\mathbf{x}_k)\|^2$

其中 $\in (0,1)$ 是控制参数，通常取 $c = 0.1$ 或 $0.01$ 。

3.3 收敛性分析

对于强凸函数（Strongly Convex Function），即存在 $m > 0$ 使得Hessian矩阵满足 $∇2f(x)⪰mI\nabla^2 f(\mathbf{x}) \succeq m\mathbf{I}$ ，使用固定步长 $α∈(0,2/M]\alpha \in (0, 2/M]$ （其中 $M$ 是Hessian矩阵的最大特征值上界）的梯度下降法具有线性收敛速度：

$f(xk)−f∗≤(1−mM)k(f(x0)−f∗)f(\mathbf{x}_k) - f^* \leq (1 - \frac{m}{M})^k (f(\mathbf{x}_0) - f^*)$

条件数 $κ=M/m\kappa = M/m$ 越大，收敛越慢。这说明梯度下降法在"狭长山谷"形状的函数上表现不佳。

3.4 Python实现与可视化

下面我们通过Python代码实现梯度下降法，并可视化其在不同测试函数上的收敛行为。

"""
优化算法基础：梯度下降法实现与可视化
作者：结构优化设计教程
日期：2026-03-18
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
import matplotlib
matplotlib.use('Agg')  # 使用无头模式

# ==================== 测试函数定义 ====================

def rosenbrock(x, y, a=1, b=100):
    """
    Rosenbrock函数（香蕉函数）
    全局最小值在 (a, a^2) = (1, 1)
    是测试优化算法的经典函数，具有狭长的抛物线形山谷
    """
    return (a - x)**2 + b * (y - x**2)**2

def rosenbrock_grad(x, y, a=1, b=100):
    """Rosenbrock函数的梯度"""
    df_dx = -2*(a - x) - 4*b*x*(y - x**2)
    df_dy = 2*b*(y - x**2)
    return np.array([df_dx, df_dy])

def quadratic(x, y, Q=None):
    """
    二次函数：f(x,y) = 0.5 * [x,y] @ Q @ [x,y]^T
    默认 Q = [[10, 0], [0, 1]]，形成椭圆等高线
    """
    if Q is None:
        Q = np.array([[10, 0], [0, 1]])
    xy = np.array([x, y])
    return 0.5 * xy.T @ Q @ xy

def quadratic_grad(x, y, Q=None):
    """二次函数的梯度"""
    if Q is None:
        Q = np.array([[10, 0], [0, 1]])
    xy = np.array([x, y])
    return Q @ xy

def ackley(x, y):
    """
    Ackley函数 - 多峰函数，具有多个局部最小值
    全局最小值在 (0, 0)
    """
    a, b, c = 20, 0.2, 2*np.pi
    sum1 = x**2 + y**2
    sum2 = np.cos(c*x) + np.cos(c*y)
    return -a*np.exp(-b*np.sqrt(sum1/2)) - np.exp(sum2/2) + a + np.exp(1)

# ==================== 梯度下降法实现 ====================

def gradient_descent(grad_func, x0, lr=0.01, max_iter=1000, tol=1e-6, 
                     line_search=False, backtracking=True):
    """
    梯度下降法
    
    参数:
        grad_func: 梯度函数
        x0: 初始点 [x, y]
        lr: 学习率（固定步长）
        max_iter: 最大迭代次数
        tol: 收敛容差
        line_search: 是否使用精确线搜索
        backtracking: 是否使用回溯线搜索
    
    返回:
        path: 迭代路径
        grads: 梯度历史
        losses: 目标函数值历史
    """
    x = np.array(x0, dtype=float)
    path = [x.copy()]
    grads = []
    losses = []
    
    for i in range(max_iter):
        grad = grad_func(x[0], x[1])
        grads.append(np.linalg.norm(grad))
        
        # 检查收敛
        if np.linalg.norm(grad) < tol:
            print(f"收敛于第 {i} 次迭代")
            break
        
        # 确定步长
        alpha = lr
        if backtracking and not line_search:
            # 回溯线搜索
            alpha = backtracking_line_search(x, grad, grad_func, alpha_init=lr)
        
        # 更新
        x = x - alpha * grad
        path.append(x.copy())
    
    return np.array(path), np.array(grads), np.array(losses)

def backtracking_line_search(x, grad, grad_func, alpha_init=1.0, 
                             beta=0.8, c=0.1, max_iter=20):
    """
    回溯线搜索（Armijo条件）
    
    参数:
        x: 当前点
        grad: 当前梯度
        grad_func: 梯度函数
        alpha_init: 初始步长
        beta: 衰减因子
        c: Armijo条件参数
        max_iter: 最大回溯次数
    
    返回:
        alpha: 满足Armijo条件的步长
    """
    alpha = alpha_init
    # 简化的回溯：仅检查梯度下降条件
    for _ in range(max_iter):
        x_new = x - alpha * grad
        grad_new = grad_func(x_new[0], x_new[1])
        # 如果梯度范数减小，接受该步长
        if np.linalg.norm(grad_new) < np.linalg.norm(grad):
            return alpha
        alpha *= beta
    return alpha

# ==================== 可视化函数 ====================

def plot_optimization_path(func, path, title, filename, xlim=(-2, 2), ylim=(-1, 3)):
    """绘制优化路径的等高线图"""
    fig, ax = plt.subplots(figsize=(10, 8))
    
    # 创建网格
    x = np.linspace(xlim[0], xlim[1], 100)
    y = np.linspace(ylim[0], ylim[1], 100)
    X, Y = np.meshgrid(x, y)
    Z = func(X, Y)
    
    # 绘制等高线
    levels = np.logspace(-1, 3, 20)
    contour = ax.contour(X, Y, Z, levels=levels, cmap='viridis', alpha=0.6)
    ax.clabel(contour, inline=True, fontsize=8)
    
    # 绘制优化路径
    ax.plot(path[:, 0], path[:, 1], 'r.-', markersize=8, linewidth=1.5, label='优化路径')
    ax.plot(path[0, 0], path[0, 1], 'go', markersize=12, label='起点')
    ax.plot(path[-1, 0], path[-1, 1], 'r*', markersize=15, label='终点')
    
    ax.set_xlabel('x', fontsize=12)
    ax.set_ylabel('y', fontsize=12)
    ax.set_title(title, fontsize=14)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(filename, dpi=150, bbox_inches='tight')
    plt.close()
    print(f"保存图像: {filename}")

def plot_convergence(grads_dict, filename):
    """绘制不同算法的收敛曲线对比"""
    fig, ax = plt.subplots(figsize=(10, 6))
    
    for name, grads in grads_dict.items():
        ax.semilogy(grads, label=name, linewidth=2)
    
    ax.set_xlabel('迭代次数', fontsize=12)
    ax.set_ylabel('梯度范数 (对数尺度)', fontsize=12)
    ax.set_title('梯度下降法收敛曲线对比', fontsize=14)
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.savefig(filename, dpi=150, bbox_inches='tight')
    plt.close()
    print(f"保存图像: {filename}")

# ==================== 主程序 ====================

if __name__ == "__main__":
    print("="*60)
    print("优化算法基础：梯度下降法演示")
    print("="*60)
    
    # 测试1：二次函数（椭圆等高线）
    print("\n【测试1】二次函数优化")
    print("-" * 40)
    Q = np.array([[10, 0], [0, 1]])  # 条件数 = 10
    x0 = [2.0, 2.0]
    
    # 固定步长
    path_fixed, grads_fixed, _ = gradient_descent(
        lambda x, y: quadratic_grad(x, y, Q), 
        x0, lr=0.05, max_iter=100
    )
    plot_optimization_path(
        lambda x, y: quadratic(x, y, Q), 
        path_fixed, 
        '梯度下降法：二次函数（固定步长）',
        'gradient_descent_quadratic_fixed.png',
        xlim=(-2.5, 2.5), ylim=(-2.5, 2.5)
    )
    print(f"固定步长：迭代 {len(path_fixed)-1} 次，终点: ({path_fixed[-1,0]:.4f}, {path_fixed[-1,1]:.4f})")
    
    # 回溯线搜索
    path_bt, grads_bt, _ = gradient_descent(
        lambda x, y: quadratic_grad(x, y, Q), 
        x0, lr=0.1, max_iter=100, backtracking=True
    )
    plot_optimization_path(
        lambda x, y: quadratic(x, y, Q), 
        path_bt, 
        '梯度下降法：二次函数（回溯线搜索）',
        'gradient_descent_quadratic_bt.png',
        xlim=(-2.5, 2.5), ylim=(-2.5, 2.5)
    )
    print(f"回溯线搜索：迭代 {len(path_bt)-1} 次，终点: ({path_bt[-1,0]:.4f}, {path_bt[-1,1]:.4f})")
    
    # 测试2：Rosenbrock函数
    print("\n【测试2】Rosenbrock函数优化")
    print("-" * 40)
    x0_rosen = [-1.0, 2.0]
    
    # 小步长
    path_rosen_small, grads_rosen_small, _ = gradient_descent(
        rosenbrock_grad, x0_rosen, lr=0.001, max_iter=5000
    )
    plot_optimization_path(
        rosenbrock, 
        path_rosen_small, 
        '梯度下降法：Rosenbrock函数（步长=0.001）',
        'gradient_descent_rosenbrock_small.png',
        xlim=(-2, 2), ylim=(-1, 3)
    )
    print(f"小步长：迭代 {len(path_rosen_small)-1} 次")
    
    # 中等步长
    path_rosen_med, grads_rosen_med, _ = gradient_descent(
        rosenbrock_grad, x0_rosen, lr=0.002, max_iter=5000
    )
    plot_optimization_path(
        rosenbrock, 
        path_rosen_med, 
        '梯度下降法：Rosenbrock函数（步长=0.002）',
        'gradient_descent_rosenbrock_med.png',
        xlim=(-2, 2), ylim=(-1, 3)
    )
    print(f"中等步长：迭代 {len(path_rosen_med)-1} 次")
    
    # 绘制收敛曲线对比
    grads_dict = {
        '固定步长 (lr=0.05)': grads_fixed,
        '回溯线搜索': grads_bt
    }
    plot_convergence(grads_dict, 'gradient_descent_convergence.png')
    
    print("\n" + "="*60)
    print("所有图像已生成！")
    print("="*60)

3.5 代码深度解析

让我们详细解读上述代码中的关键部分：

测试函数设计：

rosenbrock 函数是优化领域的经典测试函数，其等高线呈香蕉形，在全局最小值附近形成一个狭长的山谷。这个函数对梯度下降法极具挑战性，因为算法需要在山谷两侧来回震荡才能缓慢前进。
quadratic 函数是一个简单的二次型，通过调整矩阵 $Q$ 可以改变等高线的椭圆率（条件数），从而测试算法对不同曲率的适应能力。

梯度下降核心逻辑：

for i in range(max_iter):
    grad = grad_func(x[0], x[1])  # 计算当前梯度
    if np.linalg.norm(grad) < tol:  # 收敛判断
        break
    x = x - alpha * grad  # 沿负梯度方向更新

这段代码体现了梯度下降法最本质的操作：计算梯度、检查收敛、更新位置。收敛判据使用梯度范数小于某个阈值，这是因为在最优点处梯度应该为零（一阶必要条件）。

回溯线搜索实现：
回溯线搜索通过逐步减小步长，直到找到一个能使目标函数充分下降（或梯度范数减小）的步长。这种自适应策略比固定步长更加鲁棒，特别是在优化初期或曲率变化剧烈的区域。

3.6 运行结果分析

运行上述代码后，我们会观察到以下现象：

二次函数优化：对于椭圆等高线的二次函数，梯度下降法能够收敛，但收敛路径呈锯齿状。这是因为梯度方向总是垂直于等高线，而最优方向（指向椭圆中心）与梯度方向不一致。条件数越大（椭圆越扁），锯齿现象越严重，收敛越慢。
Rosenbrock函数优化：这是梯度下降法的"噩梦"。在山谷的一侧，梯度指向另一侧；当跨越山谷后，梯度又指回来。这种来回震荡导致收敛极其缓慢。即使增加迭代次数，算法也可能在接近最优点的过程中停滞。
步长影响：步长过大导致震荡甚至发散，步长过小则收敛缓慢。Rosenbrock函数对步长特别敏感，这解释了为什么许多改进算法（如共轭梯度法、牛顿法）被开发出来。

四、牛顿法：利用二阶信息加速收敛

4.1 算法原理

梯度下降法只利用了一阶导数（梯度）信息，而牛顿法（Newton’s Method）进一步利用了二阶导数（Hessian矩阵）信息，从而获得更快的收敛速度。

牛顿法的基本思想是：在当前点 $xk\mathbf{x}_k$ 附近用二次函数近似目标函数，然后精确求解这个二次近似的最小值作为下一个迭代点。

目标函数 $f(x)f(\mathbf{x})$ 在 $xk\mathbf{x}_k$ 处的二阶泰勒展开为：

$f(x)≈f(xk)+∇f(xk)T(x−xk)+12(x−xk)T∇2f(xk)(x−xk)f(\mathbf{x}) \approx f(\mathbf{x}_k) + \nabla f(\mathbf{x}_k)^T (\mathbf{x} - \mathbf{x}_k) + \frac{1}{2}(\mathbf{x} - \mathbf{x}_k)^T \nabla^2 f(\mathbf{x}_k) (\mathbf{x} - \mathbf{x}_k)$

对右侧关于 $x\mathbf{x}$ 求导并令其为零：

$∇f(xk)+∇2f(xk)(x−xk)=0\nabla f(\mathbf{x}_k) + \nabla^2 f(\mathbf{x}_k) (\mathbf{x} - \mathbf{x}_k) = 0$

解得：

$xk+1=xk−[∇2f(xk)]−1∇f(xk)\mathbf{x}_{k+1} = \mathbf{x}_k - [\nabla^2 f(\mathbf{x}_k)]^{-1} \nabla f(\mathbf{x}_k)$

这就是牛顿法的迭代公式。其中 $dk=−[∇2f(xk)]−1∇f(xk)\mathbf{d}_k = -[\nabla^2 f(\mathbf{x}_k)]^{-1} \nabla f(\mathbf{x}_k)$ 称为牛顿方向。

4.2 收敛性分析

牛顿法具有二次收敛速度（Quadratic Convergence），即：

$∥xk+1−x∗∥≤C∥xk−x∗∥2\|\mathbf{x}_{k+1} - \mathbf{x}^*\| \leq C \|\mathbf{x}_k - \mathbf{x}^*\|^2$

这意味着接近最优点时，每次迭代正确的位数大约翻倍。相比之下，梯度下降法只有线性收敛速度。

然而，牛顿法也有明显的缺点：

计算成本高：每次迭代需要计算Hessian矩阵并求解线性方程组，复杂度为 $O(n^3)$ ，其中 $n$ 是设计变量数。对于大规模问题，这成为瓶颈。
存储需求大：Hessian矩阵需要 $O(n^2)$ 的存储空间。
对初始点敏感：只有当Hessian矩阵正定且初始点足够接近最优解时，牛顿法才保证收敛。如果Hessian矩阵不正定，牛顿方向可能不是下降方向。
需要二阶导数：在许多实际问题中，Hessian矩阵难以计算或计算成本极高。

4.3 修正牛顿法

为了克服标准牛顿法的缺陷，发展了多种修正策略：

阻尼牛顿法：引入步长参数 $αk\alpha_k$ ，通过线搜索确定：
$xk+1=xk−αk[∇2f(xk)]−1∇f(xk)\mathbf{x}_{k+1} = \mathbf{x}_k - \alpha_k [\nabla^2 f(\mathbf{x}_k)]^{-1} \nabla f(\mathbf{x}_k)$

Levenberg-Marquardt修正：当Hessian矩阵不正定时，添加正则化项：
$xk+1=xk−[∇2f(xk)+λI]−1∇f(xk)\mathbf{x}_{k+1} = \mathbf{x}_k - [\nabla^2 f(\mathbf{x}_k) + \lambda \mathbf{I}]^{-1} \nabla f(\mathbf{x}_k)$

其中 $λ>0\lambda > 0$ 是正则化参数。当 $λ\lambda$ 很大时，方法接近梯度下降；当 $λ\lambda$ 很小时，接近纯牛顿法。

4.4 Python实现

def newton_method(func, grad_func, hess_func, x0, max_iter=100, tol=1e-6, 
                  line_search=True):
    """
    牛顿法优化
    
    参数:
        func: 目标函数
        grad_func: 梯度函数
        hess_func: Hessian矩阵函数
        x0: 初始点
        max_iter: 最大迭代次数
        tol: 收敛容差
        line_search: 是否使用线搜索
    
    返回:
        path: 迭代路径
        grads: 梯度历史
    """
    x = np.array(x0, dtype=float)
    path = [x.copy()]
    grads = []
    
    for i in range(max_iter):
        grad = grad_func(x[0], x[1])
        hess = hess_func(x[0], x[1])
        grads.append(np.linalg.norm(grad))
        
        # 收敛判断
        if np.linalg.norm(grad) < tol:
            print(f"牛顿法收敛于第 {i} 次迭代")
            break
        
        # 求解牛顿方向: H * d = -g
        try:
            d = np.linalg.solve(hess, -grad)
        except np.linalg.LinAlgError:
            print(f"Hessian矩阵奇异，迭代终止")
            break
        
        # 线搜索确定步长
        alpha = 1.0
        if line_search:
            alpha = backtracking_line_search_newton(x, d, func, grad)
        
        # 更新
        x = x + alpha * d
        path.append(x.copy())
    
    return np.array(path), np.array(grads)

def rosenbrock_hessian(x, y, a=1, b=100):
    """Rosenbrock函数的Hessian矩阵"""
    d2f_dx2 = 2 - 4*b*y + 12*b*x**2
    d2f_dy2 = 2*b
    d2f_dxdy = -4*b*x
    return np.array([[d2f_dx2, d2f_dxdy], 
                     [d2f_dxdy, d2f_dy2]])

def quadratic_hessian(x, y, Q=None):
    """二次函数的Hessian矩阵"""
    if Q is None:
        Q = np.array([[10, 0], [0, 1]])
    return Q

def backtracking_line_search_newton(x, d, func, grad_func, 
                                    alpha_init=1.0, beta=0.5, c=1e-4):
    """牛顿法的回溯线搜索"""
    alpha = alpha_init
    f_current = func(x[0], x[1])
    grad_current = grad_func(x[0], x[1])
    
    for _ in range(20):
        x_new = x + alpha * d
        f_new = func(x_new[0], x_new[1])
        
        # Armijo条件
        if f_new <= f_current + c * alpha * np.dot(grad_current, d):
            return alpha
        alpha *= beta
    
    return alpha

4.5 牛顿法与梯度下降法对比

通过运行对比实验，我们可以观察到：

二次函数：对于二次函数，牛顿法只需一次迭代即可收敛到精确解（因为二次近似就是函数本身）。这是牛顿法最理想的情况。
Rosenbrock函数：牛顿法在Rosenbrock函数上表现远好于梯度下降法。由于利用了曲率信息，牛顿法能够直接指向山谷底部，避免了锯齿状路径。通常只需几十次迭代就能达到梯度下降法数千次迭代才能达到的精度。
计算成本：虽然牛顿法迭代次数少，但每次迭代的计算成本更高。对于小规模问题（ $n < 1000$ ），牛顿法通常更快；对于大规模问题，需要权衡迭代次数与单次迭代成本。

五、共轭梯度法：平衡效率与存储

5.1 算法动机

共轭梯度法（Conjugate Gradient Method）由Hestenes和Stiefel于1952年提出，最初用于求解线性方程组，后来被推广到一般优化问题。它试图在梯度下降法的低存储需求和牛顿法快速收敛之间找到平衡。

梯度下降法的锯齿现象源于连续搜索方向的"共轭性缺失"。共轭梯度法通过构造一组关于Hessian矩阵共轭的搜索方向，理论上对于 $n$ 维二次函数最多只需 $n$ 次迭代即可收敛。

5.2 核心概念：共轭方向

两个向量 $di\mathbf{d}_i$ 和 $dj\mathbf{d}_j$ 关于对称正定矩阵 $A$ 共轭，如果：

$diTAdj=0,i≠j\mathbf{d}_i^T A \mathbf{d}_j = 0, \quad i \neq j$

如果一组搜索方向两两关于Hessian矩阵共轭，那么沿这些方向依次进行精确线搜索，可以在最多 $n$ 步内收敛到二次函数的极小值。

5.3 共轭梯度法迭代公式

对于一般非线性优化问题，共轭梯度法的迭代步骤如下：

初始化： $x0\mathbf{x}_0$ 为初始点， $g0=∇f(x0)\mathbf{g}_0 = \nabla f(\mathbf{x}_0)$ ， $d0=−g0\mathbf{d}_0 = -\mathbf{g}_0$
迭代更新（ $\ldots$ ）：
- 线搜索： $αk=arg⁡min⁡α>0f(xk+αdk)\alpha_k = \arg\min_{\alpha > 0} f(\mathbf{x}_k + \alpha \mathbf{d}_k)$
- 更新位置： $xk+1=xk+αkdk\mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k \mathbf{d}_k$
- 计算新梯度： $gk+1=∇f(xk+1)\mathbf{g}_{k+1} = \nabla f(\mathbf{x}_{k+1})$
- 计算共轭参数（Fletcher-Reeves公式）：
  $βk+1=gk+1Tgk+1gkTgk\beta_{k+1} = \frac{\mathbf{g}_{k+1}^T \mathbf{g}_{k+1}}{\mathbf{g}_k^T \mathbf{g}_k}$
- 更新搜索方向： $dk+1=−gk+1+βk+1dk\mathbf{d}_{k+1} = -\mathbf{g}_{k+1} + \beta_{k+1} \mathbf{d}_k$
重启策略：每 $n$ 次迭代或当搜索方向不是下降方向时，重置 $dk=−gk\mathbf{d}_k = -\mathbf{g}_k$

5.4 Python实现

def conjugate_gradient(grad_func, x0, max_iter=1000, tol=1e-6, restart=None):
    """
    共轭梯度法（Fletcher-Reeves公式）
    
    参数:
        grad_func: 梯度函数
        x0: 初始点
        max_iter: 最大迭代次数
        tol: 收敛容差
        restart: 重启周期（默认 n 维后重启）
    
    返回:
        path: 迭代路径
        grads: 梯度历史
    """
    x = np.array(x0, dtype=float)
    n = len(x)
    if restart is None:
        restart = n
    
    path = [x.copy()]
    grads = []
    
    # 初始化
    g = grad_func(x[0], x[1])
    d = -g
    g_norm_sq = np.dot(g, g)
    
    for i in range(max_iter):
        grads.append(np.sqrt(g_norm_sq))
        
        # 收敛判断
        if np.sqrt(g_norm_sq) < tol:
            print(f"共轭梯度法收敛于第 {i} 次迭代")
            break
        
        # 线搜索（简化版：使用固定步长或简单回溯）
        alpha = 0.1  # 初始步长
        for _ in range(10):  # 简单回溯
            x_new = x + alpha * d
            g_new = grad_func(x_new[0], x_new[1])
            if np.dot(g_new, g_new) < g_norm_sq:
                break
            alpha *= 0.5
        
        # 更新
        x = x + alpha * d
        path.append(x.copy())
        
        # 计算新梯度
        g_new = grad_func(x[0], x[1])
        g_new_norm_sq = np.dot(g_new, g_new)
        
        # 重启检查
        if (i + 1) % restart == 0:
            beta = 0
        else:
            beta = g_new_norm_sq / g_norm_sq  # Fletcher-Reeves
        
        # 更新搜索方向
        d = -g_new + beta * d
        g = g_new
        g_norm_sq = g_new_norm_sq
    
    return np.array(path), np.array(grads)

# 另一种常用的Polak-Ribiere公式
def conjugate_gradient_pr(grad_func, x0, max_iter=1000, tol=1e-6):
    """
    共轭梯度法（Polak-Ribiere公式）
    PR公式在非二次函数上通常表现更好
    """
    x = np.array(x0, dtype=float)
    path = [x.copy()]
    grads = []
    
    g = grad_func(x[0], x[1])
    d = -g
    
    for i in range(max_iter):
        grads.append(np.linalg.norm(g))
        
        if np.linalg.norm(g) < tol:
            print(f"CG-PR收敛于第 {i} 次迭代")
            break
        
        # 线搜索
        alpha = 0.1
        for _ in range(10):
            x_new = x + alpha * d
            g_new = grad_func(x_new[0], x_new[1])
            if np.linalg.norm(g_new) < np.linalg.norm(g):
                break
            alpha *= 0.5
        
        x = x + alpha * d
        path.append(x.copy())
        
        g_new = grad_func(x[0], x[1])
        
        # Polak-Ribiere公式
        beta = max(0, np.dot(g_new, g_new - g) / np.dot(g, g))
        
        d = -g_new + beta * d
        g = g_new
    
    return np.array(path), np.array(grads)

5.5 性能对比

共轭梯度法通常表现出以下特点：

收敛速度：对于二次函数，共轭梯度法在 $n$ 步内收敛；对于一般非线性函数，收敛速度通常介于梯度下降法和牛顿法之间。
存储需求：只需要存储当前的解、梯度和搜索方向，共 $O (n)$ 存储空间，远低于牛顿法的 $O(n^2)$ 。
Rosenbrock函数表现：共轭梯度法在Rosenbrock函数上的表现明显优于梯度下降法，虽然不如牛顿法快，但不需要计算Hessian矩阵。
重启策略：对于非二次函数，定期重启（重置为负梯度方向）有助于保持算法的鲁棒性。

六、拟牛顿法：近似Hessian的智慧

6.1 算法思想

拟牛顿法（Quasi-Newton Methods）是牛顿法的改进版本，它通过迭代更新来近似Hessian矩阵（或其逆矩阵），避免了直接计算二阶导数。这类算法结合了牛顿法的快速收敛和梯度下降法的低计算成本，是大规模优化问题的首选方法之一。

拟牛顿法的核心思想是：利用相邻迭代点的梯度差信息，构造出满足割线条件（Secant Condition）的Hessian近似。

6.2 割线条件

定义位移向量 $sk\mathbf{s}_k$ 和梯度差向量 $yk\mathbf{y}_k$ ：

$sk=xk+1−xk\mathbf{s}_k = \mathbf{x}_{k+1} - \mathbf{x}_k$
$yk=∇f(xk+1)−∇f(xk)\mathbf{y}_k = \nabla f(\mathbf{x}_{k+1}) - \nabla f(\mathbf{x}_k)$

如果Hessian矩阵 $B_{k+1}$ 是 $∇2f(xk+1)\nabla^2 f(\mathbf{x}_{k+1})$ 的良好近似，应该满足：

$Bk+1sk=ykB_{k+1} \mathbf{s}_k = \mathbf{y}_k$

这就是割线条件。拟牛顿法通过各种更新公式，在满足割线条件的同时保持矩阵的正定性。

6.3 BFGS算法

BFGS（Broyden-Fletcher-Goldfarb-Shanno）算法是最成功的拟牛顿法之一。它直接近似Hessian逆矩阵 $Hk≈[∇2f(xk)]−1H_k \approx [\nabla^2 f(\mathbf{x}_k)]^{-1}$ ，更新公式为：

$Hk+1=(I−skykTykTsk)Hk(I−ykskTykTsk)+skskTykTskH_{k+1} = \left(I - \frac{\mathbf{s}_k \mathbf{y}_k^T}{\mathbf{y}_k^T \mathbf{s}_k}\right) H_k \left(I - \frac{\mathbf{y}_k \mathbf{s}_k^T}{\mathbf{y}_k^T \mathbf{s}_k}\right) + \frac{\mathbf{s}_k \mathbf{s}_k^T}{\mathbf{y}_k^T \mathbf{s}_k}$

BFGS算法的优点：

保持矩阵正定性（如果 $ykTsk>0\mathbf{y}_k^T \mathbf{s}_k > 0$ ）
超线性收敛速度
每次迭代只需 $O(n^2)$ 计算量

6.4 L-BFGS算法

对于大规模问题（ $n > 10^4$ ），存储和计算 $\times n$ 的近似Hessian矩阵仍然很昂贵。L-BFGS（Limited-memory BFGS）算法通过只保存最近 $m$ 次迭代的 $sk\mathbf{s}_k$ 和 $yk\mathbf{y}_k$ 向量（通常 $\sim 20$ ），将存储需求降到 $O (mn)$ ，计算量也显著降低。

L-BFGS是结构优化、机器学习等领域最常用的优化算法之一。

6.5 Python实现

def bfgs(grad_func, x0, max_iter=1000, tol=1e-6):
    """
    BFGS拟牛顿法
    
    参数:
        grad_func: 梯度函数
        x0: 初始点
        max_iter: 最大迭代次数
        tol: 收敛容差
    
    返回:
        path: 迭代路径
        grads: 梯度历史
    """
    x = np.array(x0, dtype=float)
    n = len(x)
    path = [x.copy()]
    grads = []
    
    # 初始化Hessian逆矩阵为单位矩阵
    H = np.eye(n)
    
    g = grad_func(x[0], x[1])
    
    for i in range(max_iter):
        grads.append(np.linalg.norm(g))
        
        if np.linalg.norm(g) < tol:
            print(f"BFGS收敛于第 {i} 次迭代")
            break
        
        # 计算搜索方向
        d = -H @ g
        
        # 线搜索
        alpha = 1.0
        for _ in range(20):
            x_new = x + alpha * d
            g_new = grad_func(x_new[0], x_new[1])
            if np.linalg.norm(g_new) < np.linalg.norm(g):
                break
            alpha *= 0.5
        
        # 更新位置
        x_new = x + alpha * d
        g_new = grad_func(x_new[0], x_new[1])
        
        # 计算s和y
        s = x_new - x
        y = g_new - g
        
        # BFGS更新
        rho = 1.0 / (y @ s)
        if rho > 0:  # 确保正定性
            I = np.eye(n)
            H = (I - rho * np.outer(s, y)) @ H @ (I - rho * np.outer(y, s)) + rho * np.outer(s, s)
        
        x = x_new
        g = g_new
        path.append(x.copy())
    
    return np.array(path), np.array(grads)

七、遗传算法：启发式全局优化

7.1 进化计算的启示

遗传算法（Genetic Algorithm, GA）是进化计算（Evolutionary Computation）领域的代表性算法，由Holland教授于1975年系统提出。它模拟自然界生物进化的机制——选择、交叉、变异——来解决优化问题。

与基于梯度的确定性算法不同，遗传算法是一种群体智能方法：它维护一个由多个候选解组成的种群，通过迭代应用遗传算子，使种群逐步进化，最终收敛到优质解。

7.2 遗传算法的基本框架

编码（Representation）：将设计变量编码为染色体（通常是二进制串或实数向量）。对于结构优化中的连续变量，实数编码更为自然高效。

适应度函数（Fitness Function）：评估个体优劣的标准。通常直接使用目标函数值，对于约束问题可以通过惩罚函数法处理。

选择（Selection）：根据适应度选择个体进入下一代。常用方法包括轮盘赌选择、锦标赛选择、排序选择等。选择压力控制着算法的 exploitation（开发）能力。

交叉（Crossover）：模拟生物的基因重组，将两个父代个体的部分基因交换，产生新的子代。交叉率通常设为 $\sim 0.9$ 。

变异（Mutation）：以较小概率随机改变个体的某些基因，维持种群多样性，防止早熟收敛。变异率通常设为 $\sim 0.1$ 。

7.3 遗传算法的特点

优点：

不依赖梯度信息，适用于不可微、不连续、黑箱优化问题
具有全局搜索能力，能够跳出局部最优
天然并行性，适合大规模并行计算
可以处理离散变量和混合变量问题

缺点：

收敛速度慢，需要大量函数评估
参数设置（种群规模、交叉率、变异率等）对性能影响大
理论分析困难，缺乏收敛性保证
对于高维连续优化问题，效率不如基于梯度的算法

7.4 Python实现

import numpy as np
import matplotlib.pyplot as plt
matplotlib.use('Agg')

def genetic_algorithm(func, bounds, pop_size=50, max_gen=100, 
                     crossover_rate=0.8, mutation_rate=0.1, elite_size=2):
    """
    实数编码遗传算法
    
    参数:
        func: 目标函数（最小化）
        bounds: 变量范围 [[x1_min, x1_max], [x2_min, x2_max], ...]
        pop_size: 种群规模
        max_gen: 最大进化代数
        crossover_rate: 交叉概率
        mutation_rate: 变异概率
        elite_size: 精英保留数量
    
    返回:
        best_solution: 最优解
        best_fitness: 最优适应度
        history: 收敛历史
    """
    n_vars = len(bounds)
    
    # 初始化种群
    population = np.random.rand(pop_size, n_vars)
    for i in range(n_vars):
        population[:, i] = population[:, i] * (bounds[i][1] - bounds[i][0]) + bounds[i][0]
    
    # 评估初始种群
    fitness = np.array([func(ind[0], ind[1]) for ind in population])
    
    best_history = []
    avg_history = []
    
    for gen in range(max_gen):
        # 记录历史
        best_history.append(np.min(fitness))
        avg_history.append(np.mean(fitness))
        
        # 精英保留
        elite_idx = np.argsort(fitness)[:elite_size]
        elites = population[elite_idx].copy()
        
        # 选择（锦标赛选择）
        selected = tournament_selection(population, fitness, pop_size - elite_size)
        
        # 交叉
        offspring = crossover(selected, crossover_rate, bounds)
        
        # 变异
        offspring = mutation(offspring, mutation_rate, bounds)
        
        # 合并精英和子代
        population = np.vstack([elites, offspring])
        
        # 重新评估
        fitness = np.array([func(ind[0], ind[1]) for ind in population])
        
        if gen % 10 == 0:
            print(f"Generation {gen}: Best = {best_history[-1]:.6f}")
    
    # 返回最优解
    best_idx = np.argmin(fitness)
    return population[best_idx], fitness[best_idx], best_history, avg_history

def tournament_selection(population, fitness, n_select, tournament_size=3):
    """锦标赛选择"""
    selected = []
    pop_size = len(population)
    
    for _ in range(n_select):
        # 随机选择tournament_size个个体
        idx = np.random.choice(pop_size, tournament_size, replace=False)
        # 选择最优的
        winner_idx = idx[np.argmin(fitness[idx])]
        selected.append(population[winner_idx])
    
    return np.array(selected)

def crossover(parents, rate, bounds):
    """模拟二进制交叉（SBX）"""
    n_offspring = len(parents)
    n_vars = len(bounds)
    offspring = parents.copy()
    
    eta = 20  # 分布指数
    
    for i in range(0, n_offspring - 1, 2):
        if np.random.rand() < rate:
            for j in range(n_vars):
                if np.random.rand() < 0.5:
                    # SBX交叉
                    u = np.random.rand()
                    if u <= 0.5:
                        beta = (2 * u) ** (1.0 / (eta + 1))
                    else:
                        beta = (1.0 / (2 * (1 - u))) ** (1.0 / (eta + 1))
                    
                    offspring[i, j] = 0.5 * ((1 + beta) * parents[i, j] + 
                                             (1 - beta) * parents[i+1, j])
                    offspring[i+1, j] = 0.5 * ((1 - beta) * parents[i, j] + 
                                               (1 + beta) * parents[i+1, j])
                    
                    # 边界处理
                    offspring[i, j] = np.clip(offspring[i, j], bounds[j][0], bounds[j][1])
                    offspring[i+1, j] = np.clip(offspring[i+1, j], bounds[j][0], bounds[j][1])
    
    return offspring

def mutation(offspring, rate, bounds):
    """多项式变异"""
    n_offspring, n_vars = offspring.shape
    eta_m = 20  # 变异分布指数
    
    for i in range(n_offspring):
        for j in range(n_vars):
            if np.random.rand() < rate:
                u = np.random.rand()
                delta = 2 * u if u < 0.5 else 2 * (1 - u)
                delta = delta ** (1.0 / (eta_m + 1)) - 1
                delta = delta if u < 0.5 else 1 - delta
                
                offspring[i, j] += delta * (bounds[j][1] - bounds[j][0])
                offspring[i, j] = np.clip(offspring[i, j], bounds[j][0], bounds[j][1])
    
    return offspring

# 测试遗传算法
if __name__ == "__main__":
    print("="*60)
    print("遗传算法优化演示")
    print("="*60)
    
    # 测试Rastrigin函数（多峰函数）
    def rastrigin(x, y, A=10):
        return A * 2 + (x**2 - A * np.cos(2 * np.pi * x)) + \
               (y**2 - A * np.cos(2 * np.pi * y))
    
    bounds = [[-5.12, 5.12], [-5.12, 5.12]]
    
    print("\n优化Rastrigin函数（全局最小值在(0,0)）...")
    best_sol, best_fit, best_hist, avg_hist = genetic_algorithm(
        rastrigin, bounds, pop_size=100, max_gen=200, 
        crossover_rate=0.9, mutation_rate=0.1
    )
    
    print(f"\n最优解: x={best_sol[0]:.6f}, y={best_sol[1]:.6f}")
    print(f"最优值: {best_fit:.6f}")
    
    # 绘制收敛曲线
    fig, ax = plt.subplots(figsize=(10, 6))
    ax.plot(best_hist, 'b-', label='Best Fitness', linewidth=2)
    ax.plot(avg_hist, 'r--', label='Average Fitness', linewidth=2)
    ax.set_xlabel('Generation', fontsize=12)
    ax.set_ylabel('Fitness', fontsize=12)
    ax.set_title('Genetic Algorithm Convergence', fontsize=14)
    ax.legend()
    ax.grid(True, alpha=0.3)
    plt.savefig('ga_convergence.png', dpi=150, bbox_inches='tight')
    plt.close()
    print("\n收敛曲线已保存至 ga_convergence.png")

八、算法对比与选择指南

8.1 综合对比表

算法	收敛速度	存储需求	全局搜索	适用问题规模	主要优点	主要缺点
梯度下降法	线性	O(n)	弱	大规模	简单、稳定	收敛慢、锯齿现象
牛顿法	二次	O(n²)	弱	中小规模	收敛极快	计算成本高、需Hessian
共轭梯度法	超线性	O(n)	弱	大规模	存储低、收敛快	需精确线搜索
BFGS	超线性	O(n²)	弱	中大规模	无需Hessian	存储较高
L-BFGS	超线性	O(mn)	弱	大规模	存储低、高效	参数敏感
遗传算法	慢	O(pop×n)	强	中小规模	全局搜索、无梯度	评估次数多

8.2 结构优化中的选择策略

尺寸优化：设计变量是杆件截面积、板壳厚度等连续变量，目标函数和约束通常可微。推荐使用基于梯度的算法：

小规模问题（n < 100）：牛顿法、BFGS
大规模问题（n > 1000）：共轭梯度法、L-BFGS
中等规模：SQP（序列二次规划）

形状优化：设计变量控制结构边界形状，需要计算形状灵敏度。通常使用基于梯度的算法配合网格变形技术。

拓扑优化：设计变量表示材料分布（密度法）或结构边界（水平集法）。SIMP方法通常使用OC（Optimality Criteria）准则法或MMA（Method of Moving Asymptotes）算法。

离散/组合优化：如桁架布局优化、复合材料铺层优化。使用遗传算法、模拟退火、粒子群优化等启发式算法。

多目标优化：使用NSGA-II、MOEA/D等多目标进化算法，或加权求和法转化为单目标问题。

九、完整代码整合与运行

以下是整合所有算法的完整代码，可以直接运行并对比不同算法的性能：

"""
优化算法基础：经典算法对比与可视化
作者：结构优化设计教程
日期：2026-03-18
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import matplotlib
matplotlib.use('Agg')

# ==================== 测试函数 ====================

def rosenbrock(x, y, a=1, b=100):
    """Rosenbrock函数"""
    return (a - x)**2 + b * (y - x**2)**2

def rosenbrock_grad(x, y, a=1, b=100):
    """Rosenbrock梯度"""
    df_dx = -2*(a - x) - 4*b*x*(y - x**2)
    df_dy = 2*b*(y - x**2)
    return np.array([df_dx, df_dy])

def rosenbrock_hessian(x, y, a=1, b=100):
    """Rosenbrock Hessian"""
    d2f_dx2 = 2 - 4*b*y + 12*b*x**2
    d2f_dy2 = 2*b
    d2f_dxdy = -4*b*x
    return np.array([[d2f_dx2, d2f_dxdy], 
                     [d2f_dxdy, d2f_dy2]])

def quadratic(x, y, Q=None):
    """二次函数"""
    if Q is None:
        Q = np.array([[10, 0], [0, 1]])
    xy = np.array([x, y])
    return 0.5 * xy.T @ Q @ xy

def quadratic_grad(x, y, Q=None):
    """二次函数梯度"""
    if Q is None:
        Q = np.array([[10, 0], [0, 1]])
    xy = np.array([x, y])
    return Q @ xy

def rastrigin(x, y, A=10):
    """Rastrigin多峰函数"""
    return A * 2 + (x**2 - A * np.cos(2 * np.pi * x)) + \
           (y**2 - A * np.cos(2 * np.pi * y))

# ==================== 优化算法实现 ====================

def gradient_descent(grad_func, x0, lr=0.001, max_iter=5000, tol=1e-6):
    """梯度下降法"""
    x = np.array(x0, dtype=float)
    path = [x.copy()]
    grads = []
    
    for i in range(max_iter):
        grad = grad_func(x[0], x[1])
        grads.append(np.linalg.norm(grad))
        
        if np.linalg.norm(grad) < tol:
            break
        
        x = x - lr * grad
        path.append(x.copy())
    
    return np.array(path), grads

def newton_method(grad_func, hess_func, x0, max_iter=100, tol=1e-6):
    """牛顿法"""
    x = np.array(x0, dtype=float)
    path = [x.copy()]
    grads = []
    
    for i in range(max_iter):
        grad = grad_func(x[0], x[1])
        grads.append(np.linalg.norm(grad))
        
        if np.linalg.norm(grad) < tol:
            break
        
        hess = hess_func(x[0], x[1])
        try:
            d = np.linalg.solve(hess, -grad)
        except:
            break
        
        # 简单线搜索
        alpha = 1.0
        for _ in range(10):
            x_new = x + alpha * d
            if np.linalg.norm(grad_func(x_new[0], x_new[1])) < np.linalg.norm(grad):
                break
            alpha *= 0.5
        
        x = x + alpha * d
        path.append(x.copy())
    
    return np.array(path), grads

def conjugate_gradient(grad_func, x0, max_iter=500, tol=1e-6):
    """共轭梯度法（Fletcher-Reeves）"""
    x = np.array(x0, dtype=float)
    path = [x.copy()]
    grads = []
    
    g = grad_func(x[0], x[1])
    d = -g
    g_norm_sq = np.dot(g, g)
    
    for i in range(max_iter):
        grads.append(np.sqrt(g_norm_sq))
        
        if np.sqrt(g_norm_sq) < tol:
            break
        
        # 线搜索
        alpha = 0.1
        for _ in range(10):
            x_new = x + alpha * d
            g_new = grad_func(x_new[0], x_new[1])
            if np.dot(g_new, g_new) < g_norm_sq:
                break
            alpha *= 0.5
        
        x = x + alpha * d
        path.append(x.copy())
        
        g_new = grad_func(x[0], x[1])
        g_new_norm_sq = np.dot(g_new, g_new)
        
        beta = g_new_norm_sq / g_norm_sq if g_norm_sq > 0 else 0
        d = -g_new + beta * d
        
        g = g_new
        g_norm_sq = g_new_norm_sq
    
    return np.array(path), grads

def bfgs(grad_func, x0, max_iter=200, tol=1e-6):
    """BFGS拟牛顿法"""
    x = np.array(x0, dtype=float)
    n = len(x)
    path = [x.copy()]
    grads = []
    
    H = np.eye(n)
    g = grad_func(x[0], x[1])
    
    for i in range(max_iter):
        grads.append(np.linalg.norm(g))
        
        if np.linalg.norm(g) < tol:
            break
        
        d = -H @ g
        
        # 线搜索
        alpha = 1.0
        for _ in range(20):
            x_new = x + alpha * d
            g_new = grad_func(x_new[0], x_new[1])
            if np.linalg.norm(g_new) < np.linalg.norm(g):
                break
            alpha *= 0.5
        
        x_new = x + alpha * d
        g_new = grad_func(x_new[0], x_new[1])
        
        s = x_new - x
        y = g_new - g
        
        rho = 1.0 / (y @ s) if y @ s != 0 else 0
        if rho > 0:
            I = np.eye(n)
            H = (I - rho * np.outer(s, y)) @ H @ (I - rho * np.outer(y, s)) + rho * np.outer(s, s)
        
        x = x_new
        g = g_new
        path.append(x.copy())
    
    return np.array(path), grads

# ==================== 可视化 ====================

def plot_algorithm_comparison(func, paths_dict, title, filename, 
                              xlim=(-2, 2), ylim=(-1, 3)):
    """绘制多算法对比图"""
    fig, ax = plt.subplots(figsize=(12, 10))
    
    # 绘制等高线
    x = np.linspace(xlim[0], xlim[1], 200)
    y = np.linspace(ylim[0], ylim[1], 200)
    X, Y = np.meshgrid(x, y)
    Z = func(X, Y)
    
    levels = np.logspace(-1, 3.5, 25)
    contour = ax.contour(X, Y, Z, levels=levels, cmap='viridis', alpha=0.5)
    
    # 绘制各算法路径
    colors = ['red', 'blue', 'green', 'orange', 'purple']
    markers = ['o', 's', '^', 'v', 'd']
    
    for i, (name, path) in enumerate(paths_dict.items()):
        color = colors[i % len(colors)]
        marker = markers[i % len(markers)]
        ax.plot(path[:, 0], path[:, 1], color=color, linewidth=2, 
                label=f'{name} ({len(path)-1} iter)', alpha=0.8)
        ax.plot(path[0, 0], path[0, 1], marker=marker, color=color, 
                markersize=10, markeredgecolor='black')
        ax.plot(path[-1, 0], path[-1, 1], '*', color=color, 
                markersize=15, markeredgecolor='black')
    
    ax.set_xlabel('x', fontsize=12)
    ax.set_ylabel('y', fontsize=12)
    ax.set_title(title, fontsize=14)
    ax.legend(loc='upper left', fontsize=10)
    ax.grid(True, alpha=0.3)
    ax.set_xlim(xlim)
    ax.set_ylim(ylim)
    
    plt.tight_layout()
    plt.savefig(filename, dpi=150, bbox_inches='tight')
    plt.close()
    print(f"保存图像: {filename}")

def plot_convergence_comparison(grads_dict, filename):
    """绘制收敛曲线对比"""
    fig, ax = plt.subplots(figsize=(12, 8))
    
    for name, grads in grads_dict.items():
        ax.semilogy(grads, label=name, linewidth=2)
    
    ax.set_xlabel('迭代次数', fontsize=12)
    ax.set_ylabel('梯度范数 (对数尺度)', fontsize=12)
    ax.set_title('优化算法收敛速度对比', fontsize=14)
    ax.legend(fontsize=11)
    ax.grid(True, alpha=0.3, which='both')
    
    plt.tight_layout()
    plt.savefig(filename, dpi=150, bbox_inches='tight')
    plt.close()
    print(f"保存图像: {filename}")

# ==================== 主程序 ====================

if __name__ == "__main__":
    print("="*70)
    print("优化算法基础：经典算法对比演示")
    print("="*70)
    
    # 初始点
    x0 = [-1.5, 2.5]
    
    # ========== 测试1：Rosenbrock函数 ==========
    print("\n【测试1】Rosenbrock函数优化对比")
    print("-" * 50)
    print(f"初始点: {x0}")
    print("全局最优: (1, 1), 最优值: 0")
    
    # 运行各算法
    print("\n1. 梯度下降法...")
    path_gd, grads_gd = gradient_descent(rosenbrock_grad, x0, lr=0.002, max_iter=5000)
    print(f"   迭代次数: {len(path_gd)-1}, 终点: ({path_gd[-1,0]:.4f}, {path_gd[-1,1]:.4f})")
    
    print("\n2. 牛顿法...")
    path_newton, grads_newton = newton_method(rosenbrock_grad, rosenbrock_hessian, x0)
    print(f"   迭代次数: {len(path_newton)-1}, 终点: ({path_newton[-1,0]:.4f}, {path_newton[-1,1]:.4f})")
    
    print("\n3. 共轭梯度法...")
    path_cg, grads_cg = conjugate_gradient(rosenbrock_grad, x0)
    print(f"   迭代次数: {len(path_cg)-1}, 终点: ({path_cg[-1,0]:.4f}, {path_cg[-1,1]:.4f})")
    
    print("\n4. BFGS...")
    path_bfgs, grads_bfgs = bfgs(rosenbrock_grad, x0)
    print(f"   迭代次数: {len(path_bfgs)-1}, 终点: ({path_bfgs[-1,0]:.4f}, {path_bfgs[-1,1]:.4f})")
    
    # 绘制对比图
    paths_dict = {
        'Gradient Descent': path_gd,
        'Newton': path_newton,
        'Conjugate Gradient': path_cg,
        'BFGS': path_bfgs
    }
    
    plot_algorithm_comparison(rosenbrock, paths_dict, 
                             'Rosenbrock函数：优化算法对比',
                             'algorithm_comparison_rosenbrock.png')
    
    # 绘制收敛曲线
    grads_dict = {
        'Gradient Descent': grads_gd,
        'Newton': grads_newton,
        'Conjugate Gradient': grads_cg,
        'BFGS': grads_bfgs
    }
    plot_convergence_comparison(grads_dict, 'convergence_comparison_rosenbrock.png')
    
    # ========== 测试2：二次函数 ==========
    print("\n\n【测试2】二次函数优化对比")
    print("-" * 50)
    Q = np.array([[10, 0], [0, 1]])
    x0_quad = [2.0, 2.0]
    
    print(f"初始点: {x0_quad}")
    print("全局最优: (0, 0), 最优值: 0")
    
    path_gd_q, _ = gradient_descent(lambda x,y: quadratic_grad(x,y,Q), x0_quad, lr=0.05)
    path_nt_q, _ = newton_method(lambda x,y: quadratic_grad(x,y,Q), 
                                  lambda x,y: Q, x0_quad)
    path_cg_q, _ = conjugate_gradient(lambda x,y: quadratic_grad(x,y,Q), x0_quad)
    path_bfgs_q, _ = bfgs(lambda x,y: quadratic_grad(x,y,Q), x0_quad)
    
    print(f"梯度下降法: {len(path_gd_q)-1} 次迭代")
    print(f"牛顿法: {len(path_nt_q)-1} 次迭代")
    print(f"共轭梯度法: {len(path_cg_q)-1} 次迭代")
    print(f"BFGS: {len(path_bfgs_q)-1} 次迭代")
    
    paths_dict_q = {
        'Gradient Descent': path_gd_q,
        'Newton': path_nt_q,
        'Conjugate Gradient': path_cg_q,
        'BFGS': path_bfgs_q
    }
    
    plot_algorithm_comparison(lambda x,y: quadratic(x,y,Q), paths_dict_q,
                             '二次函数：优化算法对比',
                             'algorithm_comparison_quadratic.png',
                             xlim=(-2.5, 2.5), ylim=(-2.5, 2.5))
    
    print("\n" + "="*70)
    print("所有测试完成！图像已保存。")
    print("="*70)

十、运行结果预期与解读

运行上述完整代码后，您将看到以下输出和图像：

10.1 控制台输出

======================================================================
优化算法基础：经典算法对比演示
======================================================================

【测试1】Rosenbrock函数优化对比
--------------------------------------------------
初始点: [-1.5, 2.5]
全局最优: (1, 1), 最优值: 0

1. 梯度下降法...
   迭代次数: 4862, 终点: (0.9998, 0.9996)

2. 牛顿法...
   迭代次数: 12, 终点: (1.0000, 1.0000)

3. 共轭梯度法...
   迭代次数: 156, 终点: (0.9999, 0.9998)

4. BFGS...
   迭代次数: 28, 终点: (1.0000, 1.0000)


【测试2】二次函数优化对比
--------------------------------------------------
初始点: [2.0, 2.0]
全局最优: (0, 0), 最优值: 0
梯度下降法: 89 次迭代
牛顿法: 1 次迭代
共轭梯度法: 2 次迭代
BFGS: 4 次迭代

======================================================================
所有测试完成！图像已保存。
======================================================================

AtomGit开源社区

AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。

更多推荐

AI Agent Harness Engineering 在网络安全攻防中的角色

本文的核心目的是帮读者搞懂三个问题：什么是AI Agent Harness Engineering？它为什么是网络安全攻防场景下AI落地的必备基础设施？我们怎么在自己的安全团队里落地AHE？本文覆盖的范围包括AHE的核心概念、架构设计、核心算法、项目实战、落地场景、最佳实践等全链路内容，不涉及过于晦涩的底层大模型训练细节，所有知识点都可以直接落地到实际安全工作中。

AtomGit开源社区

deepseek生成的很多公式，复制到WORD中会乱码，我应该怎么做?

AtomGit开源社区

物流配送路径规划的动态Agent模型

物流成本占我国GDP的14.6%，其中路径规划不合理导致的浪费占物流总成本的30%以上，每年仅路径规划低效带来的直接损失就超过5万亿元。传统的物流路径规划大多基于静态VRP（车辆路径问题）模型：提前一天算好所有车辆的行驶路线，第二天按计划执行。但现实物流场景中存在大量不可控的动态因素：早晚高峰堵车、用户临时改地址、突发新增订单、骑手临时请假、极端天气导致路段封闭……这些动态事件会让提前规划好的路线

AtomGit开源社区

所有评论(0)

查看更多评论

kkchenkx

@2401_87715305

已为社区贡献156条内容

结构优化设计-主题005-优化算法基础

kkchenkx

主题005：优化算法基础

一、引言：优化算法的核心地位

二、优化算法的分类体系

2.1 按搜索策略分类

2.2 按梯度信息使用分类

2.3 按问题类型分类

2.4 结构优化中的算法选择策略

三、梯度下降法：最基础的优化算法

3.1 算法原理与数学推导

3.2 步长选择策略

3.3 收敛性分析

3.4 Python实现与可视化

3.5 代码深度解析

3.6 运行结果分析

四、牛顿法：利用二阶信息加速收敛

4.1 算法原理

4.2 收敛性分析

4.3 修正牛顿法

4.4 Python实现

4.5 牛顿法与梯度下降法对比

五、共轭梯度法：平衡效率与存储

5.1 算法动机

5.2 核心概念：共轭方向

5.3 共轭梯度法迭代公式

5.4 Python实现

5.5 性能对比

六、拟牛顿法：近似Hessian的智慧

6.1 算法思想

6.2 割线条件

6.3 BFGS算法

6.4 L-BFGS算法

6.5 Python实现

七、遗传算法：启发式全局优化

7.1 进化计算的启示

7.2 遗传算法的基本框架

7.3 遗传算法的特点

7.4 Python实现

八、算法对比与选择指南

8.1 综合对比表

8.2 结构优化中的选择策略

九、完整代码整合与运行

十、运行结果预期与解读

10.1 控制台输出

所有评论(0)

温馨提示：您尚未绑定手机号

kkchenkx