优化算法

SciPy 的 scipy.optimize 模块提供了丰富的数值优化算法，用于解决各种优化问题。无论是寻找函数的最小值、求解方程组，还是进行参数拟合，这个模块都能提供强大的工具。本章将详细介绍如何使用这些优化算法解决实际问题。

scipy.optimize 模块概述

scipy.optimize 模块包含以下主要功能：

标量函数优化（单变量和多变量）
约束优化和无约束优化
全局优化算法
方程求根
最小二乘拟合
线性规划
非线性规划

python

import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import warnings
warnings.filterwarnings('ignore')

# 设置绘图样式
plt.style.use('seaborn-v0_8')
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

# 查看 optimize 模块的主要功能
print("scipy.optimize 主要功能:")
functions = [attr for attr in dir(optimize) if not attr.startswith('_')]
print(f"总共 {len(functions)} 个函数和类")
print("部分函数:", functions[:15])

标量函数优化

1. 单变量函数优化

寻找最小值

python

# 定义测试函数
def f1(x):
    """简单的二次函数"""
    return x**2 + 4*x + 3

def f2(x):
    """复杂的多峰函数"""
    return x**4 - 3*x**3 + 2*x**2 + 2*x - 1

def f3(x):
    """三角函数"""
    return np.sin(x) + 0.1*x**2

# 可视化函数
x = np.linspace(-5, 5, 1000)

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

functions = [(f1, "二次函数"), (f2, "多峰函数"), (f3, "三角函数")]

for i, (func, name) in enumerate(functions):
    y = [func(xi) for xi in x]
    axes[i].plot(x, y, 'b-', linewidth=2)
    axes[i].set_title(name)
    axes[i].set_xlabel('x')
    axes[i].set_ylabel('f(x)')
    axes[i].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# 使用不同方法寻找最小值
print("单变量函数优化结果:")

for func, name in functions:
    print(f"\n{name}:")
    
    # 黄金分割法
    result_golden = optimize.minimize_scalar(func, method='golden')
    print(f"  黄金分割法: x={result_golden.x:.6f}, f(x)={result_golden.fun:.6f}")
    
    # Brent 方法
    result_brent = optimize.minimize_scalar(func, method='brent')
    print(f"  Brent 方法: x={result_brent.x:.6f}, f(x)={result_brent.fun:.6f}")
    
    # 有界优化
    result_bounded = optimize.minimize_scalar(func, bounds=(-3, 3), method='bounded')
    print(f"  有界优化: x={result_bounded.x:.6f}, f(x)={result_bounded.fun:.6f}")

寻找最大值

python

# 寻找最大值（通过最小化负函数）
def find_maximum(func, bounds=None):
    """寻找函数的最大值"""
    neg_func = lambda x: -func(x)
    
    if bounds:
        result = optimize.minimize_scalar(neg_func, bounds=bounds, method='bounded')
    else:
        result = optimize.minimize_scalar(neg_func, method='brent')
    
    return result.x, -result.fun

# 测试函数
def test_func(x):
    return -(x - 2)**2 + 5

x_max, f_max = find_maximum(test_func, bounds=(0, 4))
print(f"\n最大值优化:")
print(f"最大值点: x = {x_max:.6f}")
print(f"最大值: f(x) = {f_max:.6f}")

# 可视化
x = np.linspace(0, 4, 1000)
y = [test_func(xi) for xi in x]

plt.figure(figsize=(8, 6))
plt.plot(x, y, 'b-', linewidth=2, label='f(x) = -(x-2)² + 5')
plt.plot(x_max, f_max, 'ro', markersize=10, label=f'最大值 ({x_max:.3f}, {f_max:.3f})')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('寻找函数最大值')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

2. 多变量函数优化

无约束优化

python

# 定义多变量测试函数
def rosenbrock(x):
    """Rosenbrock 函数（香蕉函数）"""
    return 100 * (x[1] - x[0]**2)**2 + (1 - x[0])**2

def himmelblau(x):
    """Himmelblau 函数（多个全局最小值）"""
    return (x[0]**2 + x[1] - 11)**2 + (x[0] + x[1]**2 - 7)**2

def sphere(x):
    """球面函数"""
    return np.sum(x**2)

def rastrigin(x):
    """Rastrigin 函数（多峰函数）"""
    A = 10
    n = len(x)
    return A * n + np.sum(x**2 - A * np.cos(2 * np.pi * x))

# 可视化二维函数
def plot_2d_function(func, bounds, title, levels=20):
    x = np.linspace(bounds[0], bounds[1], 100)
    y = np.linspace(bounds[0], bounds[1], 100)
    X, Y = np.meshgrid(x, y)
    Z = np.zeros_like(X)
    
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            Z[i, j] = func([X[i, j], Y[i, j]])
    
    plt.figure(figsize=(10, 8))
    
    # 等高线图
    plt.subplot(2, 2, 1)
    contour = plt.contour(X, Y, Z, levels=levels)
    plt.colorbar(contour)
    plt.title(f'{title} - 等高线图')
    plt.xlabel('x₁')
    plt.ylabel('x₂')
    
    # 填充等高线图
    plt.subplot(2, 2, 2)
    contourf = plt.contourf(X, Y, Z, levels=levels, cmap='viridis')
    plt.colorbar(contourf)
    plt.title(f'{title} - 填充等高线图')
    plt.xlabel('x₁')
    plt.ylabel('x₂')
    
    # 3D 表面图
    ax = plt.subplot(2, 2, 3, projection='3d')
    surf = ax.plot_surface(X, Y, Z, cmap='viridis', alpha=0.8)
    ax.set_title(f'{title} - 3D 表面图')
    ax.set_xlabel('x₁')
    ax.set_ylabel('x₂')
    ax.set_zlabel('f(x₁, x₂)')
    
    # 3D 等高线图
    ax = plt.subplot(2, 2, 4, projection='3d')
    ax.contour(X, Y, Z, levels=levels)
    ax.set_title(f'{title} - 3D 等高线图')
    ax.set_xlabel('x₁')
    ax.set_ylabel('x₂')
    ax.set_zlabel('f(x₁, x₂)')
    
    plt.tight_layout()
    plt.show()
    
    return X, Y, Z

# 可视化 Rosenbrock 函数
X, Y, Z = plot_2d_function(rosenbrock, (-2, 2), 'Rosenbrock 函数')

# 使用不同优化方法
test_functions = [
    (rosenbrock, "Rosenbrock 函数", [0, 0]),
    (himmelblau, "Himmelblau 函数", [0, 0]),
    (sphere, "球面函数", [1, 1]),
    (rastrigin, "Rastrigin 函数", [0.5, 0.5])
]

optimization_methods = [
    'Nelder-Mead',
    'BFGS',
    'CG',
    'Powell',
    'L-BFGS-B'
]

print("\n多变量函数优化结果:")

for func, name, x0 in test_functions:
    print(f"\n{name} (初始点: {x0}):")
    
    for method in optimization_methods:
        try:
            result = optimize.minimize(func, x0, method=method)
            print(f"  {method:12s}: x={result.x}, f(x)={result.fun:.6f}, 迭代次数={result.nit if hasattr(result, 'nit') else 'N/A'}")
        except Exception as e:
            print(f"  {method:12s}: 优化失败 - {str(e)[:50]}")

梯度和海塞矩阵

python

# 定义带梯度的函数
def rosenbrock_with_grad(x):
    """Rosenbrock 函数及其梯度"""
    f = 100 * (x[1] - x[0]**2)**2 + (1 - x[0])**2
    
    # 梯度
    grad = np.zeros(2)
    grad[0] = -400 * x[0] * (x[1] - x[0]**2) - 2 * (1 - x[0])
    grad[1] = 200 * (x[1] - x[0]**2)
    
    return f, grad

def rosenbrock_hess(x):
    """Rosenbrock 函数的海塞矩阵"""
    hess = np.zeros((2, 2))
    hess[0, 0] = -400 * (x[1] - 3 * x[0]**2) + 2
    hess[0, 1] = -400 * x[0]
    hess[1, 0] = -400 * x[0]
    hess[1, 1] = 200
    
    return hess

# 比较使用和不使用梯度的优化效果
x0 = [-1, 1]

print("\n梯度优化比较:")

# 不使用梯度
result_no_grad = optimize.minimize(rosenbrock, x0, method='BFGS')
print(f"不使用梯度: x={result_no_grad.x}, f(x)={result_no_grad.fun:.8f}, 函数调用次数={result_no_grad.nfev}")

# 使用梯度
result_with_grad = optimize.minimize(rosenbrock, x0, method='BFGS', 
                                   jac=lambda x: rosenbrock_with_grad(x)[1])
print(f"使用梯度:   x={result_with_grad.x}, f(x)={result_with_grad.fun:.8f}, 函数调用次数={result_with_grad.nfev}")

# 使用海塞矩阵
result_with_hess = optimize.minimize(rosenbrock, x0, method='Newton-CG',
                                   jac=lambda x: rosenbrock_with_grad(x)[1],
                                   hess=rosenbrock_hess)
print(f"使用海塞矩阵: x={result_with_hess.x}, f(x)={result_with_hess.fun:.8f}, 函数调用次数={result_with_hess.nfev}")

# 可视化优化路径
def plot_optimization_path(func, x0, method='BFGS', bounds=(-2, 2)):
    """可视化优化路径"""
    # 记录优化路径
    path = []
    
    def callback(x):
        path.append(x.copy())
    
    # 优化
    result = optimize.minimize(func, x0, method=method, callback=callback)
    path = np.array(path)
    
    # 绘制等高线和优化路径
    x = np.linspace(bounds[0], bounds[1], 100)
    y = np.linspace(bounds[0], bounds[1], 100)
    X, Y = np.meshgrid(x, y)
    Z = np.zeros_like(X)
    
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            Z[i, j] = func([X[i, j], Y[i, j]])
    
    plt.figure(figsize=(10, 8))
    
    # 等高线图
    plt.contour(X, Y, Z, levels=20, alpha=0.6)
    plt.contourf(X, Y, Z, levels=20, alpha=0.3, cmap='viridis')
    
    # 优化路径
    if len(path) > 0:
        plt.plot(path[:, 0], path[:, 1], 'ro-', linewidth=2, markersize=6, 
                label=f'优化路径 ({method})')
        plt.plot(x0[0], x0[1], 'go', markersize=10, label='起始点')
        plt.plot(result.x[0], result.x[1], 'r*', markersize=15, label='最优点')
    
    plt.xlabel('x₁')
    plt.ylabel('x₂')
    plt.title(f'优化路径可视化 - {method}')
    plt.legend()
    plt.colorbar()
    plt.grid(True, alpha=0.3)
    plt.show()
    
    return result, path

# 可视化不同方法的优化路径
for method in ['Nelder-Mead', 'BFGS', 'CG']:
    result, path = plot_optimization_path(rosenbrock, [-1, 1], method)
    print(f"{method}: 迭代次数={len(path)}, 最终值={result.fun:.8f}")

约束优化

1. 等式约束

python

# 定义约束优化问题
# 目标函数: minimize f(x,y) = (x-1)² + (y-2)²
# 约束条件: g(x,y) = x + y - 3 = 0

def objective(x):
    """目标函数"""
    return (x[0] - 1)**2 + (x[1] - 2)**2

def constraint_eq(x):
    """等式约束"""
    return x[0] + x[1] - 3

# 定义约束
cons_eq = {'type': 'eq', 'fun': constraint_eq}

# 初始点
x0 = [0, 0]

# 求解约束优化问题
result_constrained = optimize.minimize(objective, x0, method='SLSQP', constraints=cons_eq)

print("等式约束优化:")
print(f"最优解: x = {result_constrained.x}")
print(f"最优值: f(x) = {result_constrained.fun:.6f}")
print(f"约束值: g(x) = {constraint_eq(result_constrained.x):.6f}")

# 可视化约束优化
x = np.linspace(-1, 4, 100)
y = np.linspace(-1, 4, 100)
X, Y = np.meshgrid(x, y)
Z = (X - 1)**2 + (Y - 2)**2

plt.figure(figsize=(10, 8))

# 等高线图
contour = plt.contour(X, Y, Z, levels=20)
plt.colorbar(contour)

# 约束线
x_constraint = np.linspace(-1, 4, 100)
y_constraint = 3 - x_constraint
plt.plot(x_constraint, y_constraint, 'r-', linewidth=3, label='约束: x + y = 3')

# 无约束最优点
plt.plot(1, 2, 'go', markersize=10, label='无约束最优点 (1, 2)')

# 约束最优点
plt.plot(result_constrained.x[0], result_constrained.x[1], 'ro', markersize=10, 
         label=f'约束最优点 ({result_constrained.x[0]:.2f}, {result_constrained.x[1]:.2f})')

plt.xlabel('x₁')
plt.ylabel('x₂')
plt.title('等式约束优化可视化')
plt.legend()
plt.grid(True, alpha=0.3)
plt.axis('equal')
plt.show()

2. 不等式约束

python

# 定义不等式约束优化问题
# 目标函数: minimize f(x,y) = -x - y
# 约束条件: g₁(x,y) = x² + y² - 1 ≤ 0 (圆形区域内)
#          g₂(x,y) = -x ≤ 0 (x ≥ 0)
#          g₃(x,y) = -y ≤ 0 (y ≥ 0)

def objective_ineq(x):
    """目标函数（最大化 x + y 等价于最小化 -x - y）"""
    return -x[0] - x[1]

def constraint_circle(x):
    """圆形约束"""
    return 1 - x[0]**2 - x[1]**2

def constraint_x_positive(x):
    """x ≥ 0"""
    return x[0]

def constraint_y_positive(x):
    """y ≥ 0"""
    return x[1]

# 定义约束
cons_ineq = [
    {'type': 'ineq', 'fun': constraint_circle},
    {'type': 'ineq', 'fun': constraint_x_positive},
    {'type': 'ineq', 'fun': constraint_y_positive}
]

# 初始点
x0 = [0.5, 0.5]

# 求解约束优化问题
result_ineq = optimize.minimize(objective_ineq, x0, method='SLSQP', constraints=cons_ineq)

print("\n不等式约束优化:")
print(f"最优解: x = {result_ineq.x}")
print(f"最优值: f(x) = {result_ineq.fun:.6f}")
print(f"实际最大值: {-result_ineq.fun:.6f}")
print(f"约束检查:")
print(f"  圆形约束: {constraint_circle(result_ineq.x):.6f} (≥ 0)")
print(f"  x ≥ 0: {constraint_x_positive(result_ineq.x):.6f} (≥ 0)")
print(f"  y ≥ 0: {constraint_y_positive(result_ineq.x):.6f} (≥ 0)")

# 可视化不等式约束优化
theta = np.linspace(0, 2*np.pi, 1000)
circle_x = np.cos(theta)
circle_y = np.sin(theta)

plt.figure(figsize=(10, 8))

# 可行域
theta_feasible = np.linspace(0, np.pi/2, 100)
feasible_x = np.cos(theta_feasible)
feasible_y = np.sin(theta_feasible)
plt.fill_between(feasible_x, 0, feasible_y, alpha=0.3, color='lightblue', label='可行域')

# 约束边界
plt.plot(circle_x, circle_y, 'b-', linewidth=2, label='圆形约束: x² + y² ≤ 1')
plt.axhline(y=0, color='g', linewidth=2, label='y ≥ 0')
plt.axvline(x=0, color='g', linewidth=2, label='x ≥ 0')

# 目标函数等高线
x = np.linspace(-0.2, 1.2, 100)
y = np.linspace(-0.2, 1.2, 100)
X, Y = np.meshgrid(x, y)
Z = -X - Y  # 目标函数
contour = plt.contour(X, Y, Z, levels=10, alpha=0.6, colors='red')
plt.clabel(contour, inline=True, fontsize=8)

# 最优点
plt.plot(result_ineq.x[0], result_ineq.x[1], 'ro', markersize=10, 
         label=f'最优点 ({result_ineq.x[0]:.3f}, {result_ineq.x[1]:.3f})')

plt.xlabel('x₁')
plt.ylabel('x₂')
plt.title('不等式约束优化可视化')
plt.legend()
plt.grid(True, alpha=0.3)
plt.axis('equal')
plt.xlim(-0.2, 1.2)
plt.ylim(-0.2, 1.2)
plt.show()

3. 边界约束

python

# 使用边界约束
# 目标函数: minimize f(x,y) = x² + y²
# 边界约束: -2 ≤ x ≤ 2, -1 ≤ y ≤ 3

def objective_bounds(x):
    return x[0]**2 + x[1]**2

# 定义边界
bounds = [(-2, 2), (-1, 3)]

# 初始点
x0 = [1, 2]

# 求解边界约束优化问题
result_bounds = optimize.minimize(objective_bounds, x0, method='L-BFGS-B', bounds=bounds)

print("\n边界约束优化:")
print(f"最优解: x = {result_bounds.x}")
print(f"最优值: f(x) = {result_bounds.fun:.6f}")

# 可视化边界约束
x = np.linspace(-3, 3, 100)
y = np.linspace(-2, 4, 100)
X, Y = np.meshgrid(x, y)
Z = X**2 + Y**2

plt.figure(figsize=(10, 8))

# 等高线图
contour = plt.contour(X, Y, Z, levels=20)
plt.colorbar(contour)

# 可行域
rect = plt.Rectangle((-2, -1), 4, 4, fill=False, edgecolor='red', linewidth=3, label='可行域边界')
plt.gca().add_patch(rect)

# 最优点
plt.plot(result_bounds.x[0], result_bounds.x[1], 'ro', markersize=10, 
         label=f'最优点 ({result_bounds.x[0]:.3f}, {result_bounds.x[1]:.3f})')

# 无约束最优点
plt.plot(0, 0, 'go', markersize=10, label='无约束最优点 (0, 0)')

plt.xlabel('x₁')
plt.ylabel('x₂')
plt.title('边界约束优化可视化')
plt.legend()
plt.grid(True, alpha=0.3)
plt.axis('equal')
plt.show()

全局优化

1. 差分进化算法

python

# 使用差分进化算法进行全局优化
# 测试多峰函数

def ackley(x):
    """Ackley 函数（多峰函数）"""
    a, b, c = 20, 0.2, 2*np.pi
    d = len(x)
    sum1 = np.sum(x**2)
    sum2 = np.sum(np.cos(c * x))
    return -a * np.exp(-b * np.sqrt(sum1/d)) - np.exp(sum2/d) + a + np.e

def schwefel(x):
    """Schwefel 函数"""
    return 418.9829 * len(x) - np.sum(x * np.sin(np.sqrt(np.abs(x))))

# 可视化 Ackley 函数
X, Y, Z = plot_2d_function(ackley, (-5, 5), 'Ackley 函数', levels=30)

# 全局优化方法比较
global_functions = [
    (ackley, "Ackley 函数", (-5, 5)),
    (schwefel, "Schwefel 函数", (-500, 500)),
    (rastrigin, "Rastrigin 函数", (-5, 5))
]

print("\n全局优化算法比较:")

for func, name, bounds in global_functions:
    print(f"\n{name}:")
    
    # 差分进化
    bounds_de = [bounds] * 2  # 2维问题
    result_de = optimize.differential_evolution(func, bounds_de, seed=42)
    print(f"  差分进化:     x={result_de.x}, f(x)={result_de.fun:.6f}, 函数调用次数={result_de.nfev}")
    
    # 双重退火
    result_da = optimize.dual_annealing(func, bounds_de, seed=42)
    print(f"  双重退火:     x={result_da.x}, f(x)={result_da.fun:.6f}, 函数调用次数={result_da.nfev}")
    
    # SHGO (Simplicial Homology Global Optimization)
    try:
        result_shgo = optimize.shgo(func, bounds_de)
        print(f"  SHGO:        x={result_shgo.x}, f(x)={result_shgo.fun:.6f}, 函数调用次数={result_shgo.nfev}")
    except:
        print(f"  SHGO:        优化失败")
    
    # 比较局部优化
    x0 = [np.random.uniform(bounds[0], bounds[1]) for _ in range(2)]
    result_local = optimize.minimize(func, x0, method='BFGS')
    print(f"  局部优化(BFGS): x={result_local.x}, f(x)={result_local.fun:.6f}, 函数调用次数={result_local.nfev}")

2. 模拟退火

python

# 自定义模拟退火算法
def simulated_annealing(func, bounds, n_iter=1000, temp_init=100, temp_min=1e-8, alpha=0.95):
    """
    简单的模拟退火算法实现
    """
    # 初始解
    x_current = np.array([np.random.uniform(b[0], b[1]) for b in bounds])
    f_current = func(x_current)
    
    # 最佳解
    x_best = x_current.copy()
    f_best = f_current
    
    # 记录历史
    history = {'x': [x_current.copy()], 'f': [f_current], 'temp': [temp_init]}
    
    temp = temp_init
    
    for i in range(n_iter):
        # 生成邻域解
        x_new = x_current + np.random.normal(0, 0.1, len(x_current))
        
        # 边界处理
        for j, (low, high) in enumerate(bounds):
            x_new[j] = np.clip(x_new[j], low, high)
        
        f_new = func(x_new)
        
        # 接受准则
        if f_new < f_current or np.random.random() < np.exp(-(f_new - f_current) / temp):
            x_current = x_new
            f_current = f_new
            
            # 更新最佳解
            if f_new < f_best:
                x_best = x_new.copy()
                f_best = f_new
        
        # 降温
        temp = max(temp * alpha, temp_min)
        
        # 记录历史
        if i % 10 == 0:
            history['x'].append(x_current.copy())
            history['f'].append(f_current)
            history['temp'].append(temp)
    
    return x_best, f_best, history

# 测试自定义模拟退火
x_best, f_best, history = simulated_annealing(ackley, [(-5, 5), (-5, 5)])

print(f"\n自定义模拟退火结果:")
print(f"最优解: x = {x_best}")
print(f"最优值: f(x) = {f_best:.6f}")

# 可视化模拟退火过程
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# 优化路径
history_x = np.array(history['x'])
axes[0, 0].plot(history_x[:, 0], history_x[:, 1], 'o-', alpha=0.7)
axes[0, 0].plot(x_best[0], x_best[1], 'r*', markersize=15, label='最优点')
axes[0, 0].set_xlabel('x₁')
axes[0, 0].set_ylabel('x₂')
axes[0, 0].set_title('优化路径')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 目标函数值变化
axes[0, 1].plot(history['f'])
axes[0, 1].set_xlabel('迭代次数')
axes[0, 1].set_ylabel('目标函数值')
axes[0, 1].set_title('目标函数值变化')
axes[0, 1].grid(True, alpha=0.3)

# 温度变化
axes[1, 0].plot(history['temp'])
axes[1, 0].set_xlabel('迭代次数')
axes[1, 0].set_ylabel('温度')
axes[1, 0].set_title('温度变化')
axes[1, 0].set_yscale('log')
axes[1, 0].grid(True, alpha=0.3)

# 在等高线图上显示路径
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
X, Y = np.meshgrid(x, y)
Z = np.zeros_like(X)
for i in range(X.shape[0]):
    for j in range(X.shape[1]):
        Z[i, j] = ackley([X[i, j], Y[i, j]])

axes[1, 1].contour(X, Y, Z, levels=20, alpha=0.6)
axes[1, 1].plot(history_x[:, 0], history_x[:, 1], 'ro-', alpha=0.7, markersize=3)
axes[1, 1].plot(x_best[0], x_best[1], 'r*', markersize=15, label='最优点')
axes[1, 1].set_xlabel('x₁')
axes[1, 1].set_ylabel('x₂')
axes[1, 1].set_title('等高线图上的优化路径')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

方程求根

1. 标量方程求根

python

# 定义测试方程
def equation1(x):
    """简单多项式方程: x³ - 2x - 5 = 0"""
    return x**3 - 2*x - 5

def equation2(x):
    """超越方程: e^x - 3x = 0"""
    return np.exp(x) - 3*x

def equation3(x):
    """三角方程: sin(x) - x/2 = 0"""
    return np.sin(x) - x/2

# 可视化方程
x = np.linspace(-3, 3, 1000)

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

equations = [(equation1, "x³ - 2x - 5 = 0"), 
             (equation2, "e^x - 3x = 0"), 
             (equation3, "sin(x) - x/2 = 0")]

for i, (eq, name) in enumerate(equations):
    y = [eq(xi) for xi in x]
    axes[i].plot(x, y, 'b-', linewidth=2)
    axes[i].axhline(y=0, color='r', linestyle='--', alpha=0.7)
    axes[i].set_title(name)
    axes[i].set_xlabel('x')
    axes[i].set_ylabel('f(x)')
    axes[i].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# 使用不同方法求根
root_methods = [
    ('brentq', '布伦特方法'),
    ('bisect', '二分法'),
    ('newton', '牛顿法'),
    ('secant', '割线法')
]

print("方程求根结果:")

for eq, eq_name in equations:
    print(f"\n{eq_name}:")
    
    for method, method_name in root_methods:
        try:
            if method in ['brentq', 'bisect']:
                # 需要指定区间
                if eq == equation1:
                    root = optimize.root_scalar(eq, bracket=[1, 3], method=method)
                elif eq == equation2:
                    root = optimize.root_scalar(eq, bracket=[0, 2], method=method)
                else:  # equation3
                    root = optimize.root_scalar(eq, bracket=[0, 2], method=method)
            else:
                # 需要初始猜测
                if eq == equation1:
                    root = optimize.root_scalar(eq, x0=2, method=method)
                elif eq == equation2:
                    root = optimize.root_scalar(eq, x0=1, method=method)
                else:  # equation3
                    root = optimize.root_scalar(eq, x0=1, method=method)
            
            print(f"  {method_name:12s}: x = {root.root:.6f}, f(x) = {eq(root.root):.2e}, 迭代次数 = {root.iterations}")
        except Exception as e:
            print(f"  {method_name:12s}: 求根失败 - {str(e)[:30]}")

2. 非线性方程组

python

# 定义非线性方程组
def equations_system1(vars):
    """方程组1:
    x² + y² = 4
    x - y = 1
    """
    x, y = vars
    eq1 = x**2 + y**2 - 4
    eq2 = x - y - 1
    return [eq1, eq2]

def equations_system2(vars):
    """方程组2:
    x² - y = 1
    x + y² = 2
    """
    x, y = vars
    eq1 = x**2 - y - 1
    eq2 = x + y**2 - 2
    return [eq1, eq2]

def equations_system3(vars):
    """方程组3 (3元):
    x + y + z = 6
    x² + y² + z² = 14
    x³ + y³ + z³ = 36
    """
    x, y, z = vars
    eq1 = x + y + z - 6
    eq2 = x**2 + y**2 + z**2 - 14
    eq3 = x**3 + y**3 + z**3 - 36
    return [eq1, eq2, eq3]

# 求解非线性方程组
systems = [
    (equations_system1, [1, 1], "方程组1 (2元)"),
    (equations_system2, [1, 1], "方程组2 (2元)"),
    (equations_system3, [1, 2, 3], "方程组3 (3元)")
]

print("\n非线性方程组求解:")

for system, x0, name in systems:
    print(f"\n{name}:")
    
    # 使用不同方法
    methods = ['hybr', 'lm', 'broyden1', 'anderson']
    
    for method in methods:
        try:
            sol = optimize.root(system, x0, method=method)
            residual = np.linalg.norm(system(sol.x))
            print(f"  {method:12s}: x = {sol.x}, 残差 = {residual:.2e}, 成功 = {sol.success}")
        except Exception as e:
            print(f"  {method:12s}: 求解失败 - {str(e)[:30]}")

# 可视化二元方程组的解
def plot_system_solution(system_func, solution, title):
    """可视化二元方程组的解"""
    x = np.linspace(-3, 3, 400)
    y = np.linspace(-3, 3, 400)
    X, Y = np.meshgrid(x, y)
    
    # 计算每个方程的零等高线
    Z1 = np.zeros_like(X)
    Z2 = np.zeros_like(X)
    
    for i in range(X.shape[0]):
        for j in range(X.shape[1]):
            eqs = system_func([X[i, j], Y[i, j]])
            Z1[i, j] = eqs[0]
            Z2[i, j] = eqs[1]
    
    plt.figure(figsize=(10, 8))
    
    # 绘制零等高线
    plt.contour(X, Y, Z1, levels=[0], colors='blue', linewidths=2, label='方程1')
    plt.contour(X, Y, Z2, levels=[0], colors='red', linewidths=2, label='方程2')
    
    # 标记解
    plt.plot(solution[0], solution[1], 'go', markersize=10, 
             label=f'解 ({solution[0]:.3f}, {solution[1]:.3f})')
    
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title(title)
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.axis('equal')
    plt.show()

# 可视化前两个方程组的解
for system, x0, name in systems[:2]:
    sol = optimize.root(system, x0, method='hybr')
    if sol.success:
        plot_system_solution(system, sol.x, f'{name} 的解')

最小二乘拟合

1. 线性最小二乘

python

# 生成带噪声的数据
np.random.seed(42)
n_points = 50
x_data = np.linspace(0, 10, n_points)
y_true = 2 * x_data + 3  # 真实关系: y = 2x + 3
noise = np.random.normal(0, 2, n_points)
y_data = y_true + noise

# 线性最小二乘拟合
def linear_func(x, a, b):
    return a * x + b

# 使用 curve_fit
popt_linear, pcov_linear = optimize.curve_fit(linear_func, x_data, y_data)
a_fit, b_fit = popt_linear

print("线性最小二乘拟合:")
print(f"真实参数: a = 2, b = 3")
print(f"拟合参数: a = {a_fit:.3f}, b = {b_fit:.3f}")
print(f"参数标准误差: {np.sqrt(np.diag(pcov_linear))}")

# 计算拟合优度
y_fit = linear_func(x_data, *popt_linear)
ss_res = np.sum((y_data - y_fit)**2)
ss_tot = np.sum((y_data - np.mean(y_data))**2)
r_squared = 1 - (ss_res / ss_tot)

print(f"R² = {r_squared:.4f}")

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

# 拟合结果
plt.subplot(2, 2, 1)
plt.scatter(x_data, y_data, alpha=0.6, label='观测数据')
plt.plot(x_data, y_true, 'g-', linewidth=2, label='真实关系')
plt.plot(x_data, y_fit, 'r-', linewidth=2, label=f'拟合结果 (R²={r_squared:.3f})')
plt.xlabel('x')
plt.ylabel('y')
plt.title('线性最小二乘拟合')
plt.legend()
plt.grid(True, alpha=0.3)

# 残差图
residuals = y_data - y_fit
plt.subplot(2, 2, 2)
plt.scatter(y_fit, residuals, alpha=0.6)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('拟合值')
plt.ylabel('残差')
plt.title('残差图')
plt.grid(True, alpha=0.3)

# 残差直方图
plt.subplot(2, 2, 3)
plt.hist(residuals, bins=15, alpha=0.7, edgecolor='black')
plt.xlabel('残差')
plt.ylabel('频数')
plt.title('残差分布')
plt.grid(True, alpha=0.3)

# Q-Q 图
from scipy import stats
stats.probplot(residuals, dist="norm", plot=plt.subplot(2, 2, 4))
plt.title('残差 Q-Q 图')
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

2. 非线性最小二乘

python

# 生成非线性数据
np.random.seed(42)
x_data = np.linspace(0, 4, 50)
y_true = 3 * np.exp(-2 * x_data) + 1  # 真实关系: y = 3*exp(-2x) + 1
noise = np.random.normal(0, 0.1, len(x_data))
y_data = y_true + noise

# 定义非线性函数
def exponential_func(x, a, b, c):
    return a * np.exp(-b * x) + c

def gaussian_func(x, a, b, c):
    return a * np.exp(-(x - b)**2 / (2 * c**2))

def polynomial_func(x, a, b, c, d):
    return a * x**3 + b * x**2 + c * x + d

# 非线性拟合
functions = [
    (exponential_func, [1, 1, 1], "指数函数"),
    (gaussian_func, [1, 2, 1], "高斯函数"),
    (polynomial_func, [1, 1, 1, 1], "三次多项式")
]

print("\n非线性最小二乘拟合:")

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

# 原始数据
axes[0].scatter(x_data, y_data, alpha=0.6, label='观测数据')
axes[0].plot(x_data, y_true, 'g-', linewidth=2, label='真实关系')
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].set_title('原始数据')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

for i, (func, p0, name) in enumerate(functions):
    try:
        popt, pcov = optimize.curve_fit(func, x_data, y_data, p0=p0)
        y_fit = func(x_data, *popt)
        
        # 计算拟合优度
        ss_res = np.sum((y_data - y_fit)**2)
        ss_tot = np.sum((y_data - np.mean(y_data))**2)
        r_squared = 1 - (ss_res / ss_tot)
        
        print(f"\n{name}:")
        print(f"  拟合参数: {popt}")
        print(f"  参数标准误差: {np.sqrt(np.diag(pcov))}")
        print(f"  R² = {r_squared:.4f}")
        
        # 可视化
        axes[i+1].scatter(x_data, y_data, alpha=0.6, label='观测数据')
        axes[i+1].plot(x_data, y_true, 'g-', linewidth=2, label='真实关系')
        axes[i+1].plot(x_data, y_fit, 'r-', linewidth=2, 
                      label=f'{name} (R²={r_squared:.3f})')
        axes[i+1].set_xlabel('x')
        axes[i+1].set_ylabel('y')
        axes[i+1].set_title(f'{name} 拟合')
        axes[i+1].legend()
        axes[i+1].grid(True, alpha=0.3)
        
    except Exception as e:
        print(f"\n{name}: 拟合失败 - {str(e)}")

plt.tight_layout()
plt.show()

3. 鲁棒回归

python

# 生成带异常值的数据
np.random.seed(42)
n_points = 50
x_data = np.linspace(0, 10, n_points)
y_true = 2 * x_data + 3
noise = np.random.normal(0, 1, n_points)
y_data = y_true + noise

# 添加异常值
outlier_indices = [10, 25, 40]
for idx in outlier_indices:
    y_data[idx] += np.random.choice([-8, 8])

# 定义损失函数
def huber_loss(params, x, y, delta=1.0):
    """Huber 损失函数"""
    a, b = params
    residuals = y - (a * x + b)
    
    # Huber 损失
    loss = np.where(np.abs(residuals) <= delta,
                   0.5 * residuals**2,
                   delta * (np.abs(residuals) - 0.5 * delta))
    return np.sum(loss)

def least_squares_loss(params, x, y):
    """最小二乘损失函数"""
    a, b = params
    residuals = y - (a * x + b)
    return np.sum(residuals**2)

# 比较不同的拟合方法
print("\n鲁棒回归比较:")

# 普通最小二乘
result_ols = optimize.minimize(least_squares_loss, [1, 1], args=(x_data, y_data))
a_ols, b_ols = result_ols.x
y_fit_ols = a_ols * x_data + b_ols

# Huber 回归
result_huber = optimize.minimize(huber_loss, [1, 1], args=(x_data, y_data, 1.0))
a_huber, b_huber = result_huber.x
y_fit_huber = a_huber * x_data + b_huber

# 使用 scipy 的鲁棒方法
from scipy.optimize import least_squares

def residuals_func(params, x, y):
    a, b = params
    return y - (a * x + b)

# 普通最小二乘
result_ls = least_squares(residuals_func, [1, 1], args=(x_data, y_data))
a_ls, b_ls = result_ls.x
y_fit_ls = a_ls * x_data + b_ls

# 鲁棒最小二乘 (soft_l1 损失)
result_robust = least_squares(residuals_func, [1, 1], args=(x_data, y_data), loss='soft_l1')
a_robust, b_robust = result_robust.x
y_fit_robust = a_robust * x_data + b_robust

print(f"真实参数:     a = 2.000, b = 3.000")
print(f"普通最小二乘: a = {a_ols:.3f}, b = {b_ols:.3f}")
print(f"Huber 回归:   a = {a_huber:.3f}, b = {b_huber:.3f}")
print(f"鲁棒最小二乘: a = {a_robust:.3f}, b = {b_robust:.3f}")

# 可视化比较
plt.figure(figsize=(15, 5))

# 原始数据
plt.subplot(1, 3, 1)
plt.scatter(x_data, y_data, alpha=0.6, label='观测数据')
for idx in outlier_indices:
    plt.scatter(x_data[idx], y_data[idx], color='red', s=100, marker='x', label='异常值' if idx == outlier_indices[0] else "")
plt.plot(x_data, y_true, 'g-', linewidth=2, label='真实关系')
plt.xlabel('x')
plt.ylabel('y')
plt.title('带异常值的数据')
plt.legend()
plt.grid(True, alpha=0.3)

# 拟合结果比较
plt.subplot(1, 3, 2)
plt.scatter(x_data, y_data, alpha=0.6, label='观测数据')
for idx in outlier_indices:
    plt.scatter(x_data[idx], y_data[idx], color='red', s=100, marker='x')
plt.plot(x_data, y_true, 'g-', linewidth=2, label='真实关系')
plt.plot(x_data, y_fit_ols, 'b--', linewidth=2, label='普通最小二乘')
plt.plot(x_data, y_fit_robust, 'r-', linewidth=2, label='鲁棒回归')
plt.xlabel('x')
plt.ylabel('y')
plt.title('拟合结果比较')
plt.legend()
plt.grid(True, alpha=0.3)

# 残差比较
residuals_ols = y_data - y_fit_ols
residuals_robust = y_data - y_fit_robust

plt.subplot(1, 3, 3)
plt.scatter(y_fit_ols, residuals_ols, alpha=0.6, label='普通最小二乘')
plt.scatter(y_fit_robust, residuals_robust, alpha=0.6, label='鲁棒回归')
plt.axhline(y=0, color='k', linestyle='--')
plt.xlabel('拟合值')
plt.ylabel('残差')
plt.title('残差比较')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

线性规划

python

# 线性规划问题
# 目标函数: maximize 3x₁ + 2x₂
# 约束条件: x₁ + x₂ ≤ 4
#          2x₁ + x₂ ≤ 6
#          x₁, x₂ ≥ 0

from scipy.optimize import linprog

# 定义线性规划问题（注意：linprog 求最小值，所以目标函数系数要取负）
c = [-3, -2]  # 目标函数系数（取负号因为要最大化）
A_ub = [[1, 1],   # 不等式约束矩阵
        [2, 1]]
b_ub = [4, 6]     # 不等式约束右端
bounds = [(0, None), (0, None)]  # 变量边界

# 求解线性规划
result_lp = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')

print("\n线性规划求解结果:")
print(f"最优解: x₁ = {result_lp.x[0]:.3f}, x₂ = {result_lp.x[1]:.3f}")
print(f"最优值: {-result_lp.fun:.3f}")  # 取负号因为原问题是最大化
print(f"求解状态: {result_lp.message}")

# 可视化线性规划
x1 = np.linspace(0, 5, 400)
x2 = np.linspace(0, 5, 400)
X1, X2 = np.meshgrid(x1, x2)

plt.figure(figsize=(10, 8))

# 可行域
# 约束1: x₁ + x₂ ≤ 4
x2_constraint1 = 4 - x1
# 约束2: 2x₁ + x₂ ≤ 6
x2_constraint2 = 6 - 2*x1

# 绘制约束线
plt.plot(x1, x2_constraint1, 'b-', linewidth=2, label='x₁ + x₂ ≤ 4')
plt.plot(x1, x2_constraint2, 'r-', linewidth=2, label='2x₁ + x₂ ≤ 6')
plt.axhline(y=0, color='k', linewidth=1)
plt.axvline(x=0, color='k', linewidth=1)

# 填充可行域
feasible_x1 = np.linspace(0, 3, 100)
feasible_x2_upper = np.minimum(4 - feasible_x1, 6 - 2*feasible_x1)
feasible_x2_upper = np.maximum(feasible_x2_upper, 0)
plt.fill_between(feasible_x1, 0, feasible_x2_upper, alpha=0.3, color='lightblue', label='可行域')

# 目标函数等高线
Z = 3*X1 + 2*X2
contour = plt.contour(X1, X2, Z, levels=10, alpha=0.6, colors='green')
plt.clabel(contour, inline=True, fontsize=8)

# 最优点
plt.plot(result_lp.x[0], result_lp.x[1], 'ro', markersize=10, 
         label=f'最优点 ({result_lp.x[0]:.2f}, {result_lp.x[1]:.2f})')

plt.xlabel('x₁')
plt.ylabel('x₂')
plt.title('线性规划可视化')
plt.legend()
plt.grid(True, alpha=0.3)
plt.xlim(0, 5)
plt.ylim(0, 5)
plt.show()

实际应用案例

1. 投资组合优化

python

# 投资组合优化问题
# 假设有3只股票，历史收益率数据
np.random.seed(42)
n_assets = 3
n_periods = 252  # 一年的交易日

# 生成模拟的收益率数据
returns = np.random.multivariate_normal(
    mean=[0.08, 0.12, 0.15],  # 年化收益率
    cov=[[0.04, 0.01, 0.02],  # 协方差矩阵
         [0.01, 0.09, 0.03],
         [0.02, 0.03, 0.16]],
    size=n_periods
) / np.sqrt(252)  # 转换为日收益率

# 计算统计量
mean_returns = np.mean(returns, axis=0)
cov_matrix = np.cov(returns.T)

print("投资组合优化:")
print(f"预期收益率: {mean_returns * 252:.3f}")  # 年化
print(f"协方差矩阵:\n{cov_matrix * 252:.4f}")  # 年化

# 定义优化问题
def portfolio_variance(weights, cov_matrix):
    """投资组合方差"""
    return np.dot(weights.T, np.dot(cov_matrix, weights))

def portfolio_return(weights, mean_returns):
    """投资组合收益率"""
    return np.dot(weights, mean_returns)

# 约束条件
constraints = [
    {'type': 'eq', 'fun': lambda x: np.sum(x) - 1}  # 权重和为1
]

# 边界条件（权重在0-1之间）
bounds = tuple((0, 1) for _ in range(n_assets))

# 初始猜测
x0 = np.array([1/n_assets] * n_assets)

# 最小方差组合
result_min_var = optimize.minimize(
    portfolio_variance, x0, args=(cov_matrix,),
    method='SLSQP', bounds=bounds, constraints=constraints
)

min_var_weights = result_min_var.x
min_var_return = portfolio_return(min_var_weights, mean_returns) * 252
min_var_risk = np.sqrt(portfolio_variance(min_var_weights, cov_matrix) * 252)

print(f"\n最小方差组合:")
print(f"权重: {min_var_weights}")
print(f"年化收益率: {min_var_return:.3f}")
print(f"年化波动率: {min_var_risk:.3f}")

# 有效前沿
target_returns = np.linspace(min_var_return, 0.20, 50)
efficient_portfolios = []

for target in target_returns:
    # 添加收益率约束
    cons = [
        {'type': 'eq', 'fun': lambda x: np.sum(x) - 1},
        {'type': 'eq', 'fun': lambda x, target=target: portfolio_return(x, mean_returns) * 252 - target}
    ]
    
    result = optimize.minimize(
        portfolio_variance, x0, args=(cov_matrix,),
        method='SLSQP', bounds=bounds, constraints=cons
    )
    
    if result.success:
        risk = np.sqrt(portfolio_variance(result.x, cov_matrix) * 252)
        efficient_portfolios.append((target, risk, result.x))

# 可视化有效前沿
if efficient_portfolios:
    returns_eff = [p[0] for p in efficient_portfolios]
    risks_eff = [p[1] for p in efficient_portfolios]
    
    plt.figure(figsize=(10, 8))
    plt.plot(risks_eff, returns_eff, 'b-', linewidth=2, label='有效前沿')
    plt.scatter(min_var_risk, min_var_return, color='red', s=100, 
                label=f'最小方差组合\n(收益率={min_var_return:.3f}, 风险={min_var_risk:.3f})')
    
    # 单个资产
    for i in range(n_assets):
        asset_return = mean_returns[i] * 252
        asset_risk = np.sqrt(cov_matrix[i, i] * 252)
        plt.scatter(asset_risk, asset_return, s=100, label=f'资产 {i+1}')
    
    plt.xlabel('风险 (年化波动率)')
    plt.ylabel('收益率 (年化)')
    plt.title('投资组合有效前沿')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()

2. 生产计划优化

python

# 生产计划优化问题
# 某工厂生产两种产品，需要优化生产计划

# 问题参数
profit = [40, 30]  # 每单位产品的利润
labor_hours = [2, 1]  # 每单位产品需要的工时
material_cost = [10, 8]  # 每单位产品的材料成本

# 约束条件
max_labor = 100  # 最大工时
max_material_budget = 800  # 最大材料预算
min_product1 = 10  # 产品1的最小产量
max_product2 = 50  # 产品2的最大产量

# 定义优化问题
def production_profit(x):
    """负利润（因为要最大化利润）"""
    return -(profit[0] * x[0] + profit[1] * x[1])

# 约束函数
def labor_constraint(x):
    return max_labor - (labor_hours[0] * x[0] + labor_hours[1] * x[1])

def material_constraint(x):
    return max_material_budget - (material_cost[0] * x[0] + material_cost[1] * x[1])

def min_product1_constraint(x):
    return x[0] - min_product1

def max_product2_constraint(x):
    return max_product2 - x[1]

# 定义约束
constraints_prod = [
    {'type': 'ineq', 'fun': labor_constraint},
    {'type': 'ineq', 'fun': material_constraint},
    {'type': 'ineq', 'fun': min_product1_constraint},
    {'type': 'ineq', 'fun': max_product2_constraint}
]

# 边界条件
bounds_prod = [(0, None), (0, None)]

# 初始猜测
x0_prod = [20, 30]

# 求解
result_prod = optimize.minimize(
    production_profit, x0_prod,
    method='SLSQP', bounds=bounds_prod, constraints=constraints_prod
)

optimal_production = result_prod.x
max_profit = -result_prod.fun

print("\n生产计划优化结果:")
print(f"最优生产计划: 产品1 = {optimal_production[0]:.1f} 单位, 产品2 = {optimal_production[1]:.1f} 单位")
print(f"最大利润: {max_profit:.2f}")

# 检查约束
print(f"\n约束检查:")
print(f"工时使用: {labor_hours[0] * optimal_production[0] + labor_hours[1] * optimal_production[1]:.1f} / {max_labor}")
print(f"材料成本: {material_cost[0] * optimal_production[0] + material_cost[1] * optimal_production[1]:.1f} / {max_material_budget}")
print(f"产品1产量: {optimal_production[0]:.1f} (≥ {min_product1})")
print(f"产品2产量: {optimal_production[1]:.1f} (≤ {max_product2})")

性能优化技巧

1. 选择合适的算法

python

# 比较不同算法的性能
import time

def benchmark_optimization_methods():
    """基准测试不同优化方法的性能"""
    
    # 测试函数
    test_functions = [
        (rosenbrock, [-1, 1], "Rosenbrock"),
        (himmelblau, [0, 0], "Himmelblau"),
        (ackley, [1, 1], "Ackley")
    ]
    
    methods = ['Nelder-Mead', 'BFGS', 'CG', 'Powell', 'L-BFGS-B']
    
    results = {}
    
    for func, x0, func_name in test_functions:
        print(f"\n{func_name} 函数性能测试:")
        results[func_name] = {}
        
        for method in methods:
            try:
                start_time = time.time()
                result = optimize.minimize(func, x0, method=method)
                end_time = time.time()
                
                results[func_name][method] = {
                    'time': end_time - start_time,
                    'fun_value': result.fun,
                    'nfev': result.nfev if hasattr(result, 'nfev') else 'N/A',
                    'success': result.success
                }
                
                print(f"  {method:12s}: 时间={end_time - start_time:.4f}s, "
                      f"函数值={result.fun:.6f}, 调用次数={result.nfev if hasattr(result, 'nfev') else 'N/A'}")
                      
            except Exception as e:
                print(f"  {method:12s}: 失败 - {str(e)[:30]}")
    
    return results

# 运行基准测试
benchmark_results = benchmark_optimization_methods()

2. 并行优化

python

# 并行全局优化
from concurrent.futures import ProcessPoolExecutor
import multiprocessing

def parallel_optimization(func, bounds, n_trials=10, n_processes=None):
    """并行运行多次优化"""
    
    if n_processes is None:
        n_processes = multiprocessing.cpu_count()
    
    def single_optimization(seed):
        np.random.seed(seed)
        x0 = [np.random.uniform(b[0], b[1]) for b in bounds]
        result = optimize.minimize(func, x0, method='L-BFGS-B', bounds=bounds)
        return result.x, result.fun, result.success
    
    # 并行执行
    with ProcessPoolExecutor(max_workers=n_processes) as executor:
        futures = [executor.submit(single_optimization, i) for i in range(n_trials)]
        results = [future.result() for future in futures]
    
    # 找到最佳结果
    successful_results = [r for r in results if r[2]]  # 只考虑成功的结果
    
    if successful_results:
        best_result = min(successful_results, key=lambda x: x[1])
        return best_result[0], best_result[1], len(successful_results)
    else:
        return None, None, 0

# 测试并行优化
print("\n并行优化测试:")
start_time = time.time()
best_x, best_f, n_success = parallel_optimization(ackley, [(-5, 5), (-5, 5)], n_trials=20)
end_time = time.time()

if best_x is not None:
    print(f"最佳解: x = {best_x}")
    print(f"最佳值: f(x) = {best_f:.6f}")
    print(f"成功次数: {n_success}/20")
    print(f"总时间: {end_time - start_time:.2f}s")
else:
    print("所有优化尝试都失败了")

小结

本章详细介绍了 SciPy 的优化算法模块 scipy.optimize，包括：

标量函数优化：单变量和多变量函数的最小值/最大值求解
约束优化：等式约束、不等式约束和边界约束的处理
全局优化：差分进化、模拟退火等全局搜索算法
方程求根：标量方程和非线性方程组的求解
最小二乘拟合：线性和非线性参数估计，鲁棒回归
线性规划：线性目标函数和约束的优化问题
实际应用：投资组合优化、生产计划等实际案例
性能优化：算法选择和并行计算技巧

优化算法是科学计算和工程应用中的核心工具，掌握这些方法对于解决实际问题具有重要意义。

练习题

实现一个遗传算法来优化 Rastrigin 函数，并与 SciPy 的差分进化算法比较性能。
设计一个多目标优化问题，同时最小化成本和最大化质量，使用约束优化方法求解。
对一组实验数据进行非线性拟合，比较不同损失函数（最小二乘、Huber、绝对值）的效果。
实现一个简单的投资组合优化器，考虑交易成本和最小持仓限制。
使用线性规划解决一个运输问题：从多个供应点向多个需求点运输货物，最小化总运输成本。

优化算法 ​

scipy.optimize 模块概述 ​

标量函数优化 ​

1. 单变量函数优化 ​

寻找最小值 ​

寻找最大值 ​

2. 多变量函数优化 ​

无约束优化 ​

梯度和海塞矩阵 ​

约束优化 ​

1. 等式约束 ​

2. 不等式约束 ​

3. 边界约束 ​

全局优化 ​

1. 差分进化算法 ​

2. 模拟退火 ​

方程求根 ​

1. 标量方程求根 ​

2. 非线性方程组 ​

最小二乘拟合 ​

1. 线性最小二乘 ​

2. 非线性最小二乘 ​

3. 鲁棒回归 ​

线性规划 ​

实际应用案例 ​

1. 投资组合优化 ​

2. 生产计划优化 ​

性能优化技巧 ​

1. 选择合适的算法 ​

2. 并行优化 ​

小结 ​

练习题 ​

优化算法

scipy.optimize 模块概述

标量函数优化

1. 单变量函数优化

寻找最小值

寻找最大值

2. 多变量函数优化

无约束优化

梯度和海塞矩阵

约束优化

1. 等式约束

2. 不等式约束

3. 边界约束

全局优化

1. 差分进化算法

2. 模拟退火

方程求根

1. 标量方程求根

2. 非线性方程组

最小二乘拟合

1. 线性最小二乘

2. 非线性最小二乘

3. 鲁棒回归

线性规划

实际应用案例

1. 投资组合优化

2. 生产计划优化

性能优化技巧

1. 选择合适的算法

2. 并行优化

小结

练习题