使用 NumPy 和 SciPy 解决带线性约束的线性方程组-Python教程-PHP中文网

使用 numpy 和 scipy 解决带线性约束的线性方程组

本文探讨了如何在存在线性约束的情况下，有效求解线性方程组 AX=b。通过对比基于优化的 scipy.optimize.minimize 方法与直接的 np.linalg.lstsq 最小二乘法，阐明了将线性约束整合到方程组中并使用最小二乘求解器是处理此类问题的更优选择，尤其适用于寻求精确或最佳拟合解的场景。

引言：带约束的线性方程组求解挑战

在科学和工程计算中，我们经常需要解决形如 AX=b 的线性方程组，其中 A 是系数矩阵，X 是未知数向量，b 是常数向量。然而，实际问题往往伴随着对 X 中元素的一些额外限制，即约束条件。当这些约束是线性的时候，如何有效地将它们融入到求解过程中，并找到一个既满足原始方程组又符合所有约束的解，是一个常见的挑战。

一个常见的误区是尝试将约束条件作为惩罚项或通过优化方法来解决。虽然优化方法（如 scipy.optimize.minimize）可以处理约束，但如果目标是精确求解 AX=b 并在满足约束的同时最小化残差，那么直接的最小二乘法结合系统增广往往是更简洁和高效的方案。

传统优化方法的局限性

考虑一个典型的场景：我们有一个 8x8 的矩阵 A 和 8x1 的向量 b，需要求解 X。同时，X 的元素之间存在以下线性约束：

0.5 * (y1 + y2) = 0
0.5 * (x3 + x4) = 0
0.5 * (y3 + y4) = 0

其中 X = [x1, y1, x2, y2, x3, y3, x4, y4]。

如果尝试使用 scipy.optimize.minimize，我们通常会定义一个目标函数来最小化 ||AX - b||^2（即残差的平方和），并将约束作为等式约束传递给优化器。

import numpy as np
from scipy import optimize

# 示例数据
A = np.array([
    [-261.60,  11.26,      0.0,    0.0,     0.0,    0.0,     0.0,    0.0],
    [   4.07, -12.75,      0.0,    0.0,     0.0,    0.0,     0.0,    0.0],
    [    0.0,    0.0,  -158.63,  -5.65,     0.0,    0.0,     0.0,    0.0],
    [    0.0,    0.0,    -2.81, -12.14,     0.0,    0.0,     0.0,    0.0],
    [    0.0,    0.0,      0.0,    0.0, -265.99,  19.29,     0.0,    0.0],
    [    0.0,    0.0,      0.0,    0.0,   12.59, -12.34,     0.0,    0.0],
    [    0.0,    0.0,      0.0,    0.0,     0.0,    0.0, -166.25, -12.63],
    [    0.0,    0.0,      0.0,    0.0,     0.0,    0.0,   -8.40, -11.14]
])

b = np.array([
     -6.95,
     16.35,
     -0.96,
     16.35,
     19.19,
    -15.85,
    -12.36,
    -15.63]).reshape(-1, 1)

def objective_function(x):
    """目标函数：最小化 (AX - b) 的L2范数平方"""
    return np.sum((np.dot(A, x) - b.flatten())**2)

def constraints(x):
    """线性等式约束函数"""
    # X = [x1, y1, x2, y2, x3, y3, x4, y4]
    # 索引: x[0]=x1, x[1]=y1, x[2]=x2, x[3]=y2, x[4]=x3, x[5]=y3, x[6]=x4, x[7]=y4
    return np.array([
        0.5 * (x[1] + x[3]),  # 0.5*(y1 + y2) = 0
        0.5 * (x[4] + x[6]),  # 0.5*(x3 + x4) = 0
        0.5 * (x[5] + x[7])   # 0.5*(y3 + y4) = 0
    ])

cons = {'type': 'eq', 'fun': constraints}
res = optimize.minimize(objective_function, np.zeros(A.shape[1]), method='SLSQP', constraints=cons)

x_optimized = res.x

print("优化器找到的解 X:")
print(x_optimized)
print("\n验证约束条件 (应接近于0):")
print(constraints(x_optimized))
print("\n验证 AX 与 b 的匹配程度:")
print(np.matmul(A, x_optimized).reshape(-1, 1))
print("\n期望的 b 向量:")
print(b)

登录后复制

运行上述代码，会发现优化器虽然成功地使约束条件接近于零，但 np.matmul(A, x_optimized) 的结果与原始 b 向量仍存在显著差异。这是因为 minimize 函数的目标是找到一个 X，使得目标函数（即 ||AX - b||^2）在满足约束的前提下达到最小值。如果原始系统 AX=b 在满足约束的情况下本身就没有精确解（即 ||AX - b||^2 的最小值为非零），那么 minimize 就不会强制 AX 等于 b，而是找到一个残差最小的解。

增广系统与最小二乘法：更直接的解决方案

对于线性约束，存在一种更直接、更符合数学原理的方法：将约束条件直接整合到原始的线性方程组中，形成一个增广系统，然后使用最小二乘法求解。

每个线性约束 c_1 x_1 + c_2 x_2 + ... + c_n x_n = d 都可以被视为原始系统的一个额外方程。通过将这些约束方程添加到 A 矩阵和 b 向量中，我们构建了一个新的、可能过定（方程数多于未知数）的线性系统 A_aug X = b_aug。然后，我们可以使用 np.linalg.lstsq 来求解这个增广系统，它会找到一个 X，使得 ||A_aug X - b_aug||^2 最小。如果增广系统有精确解，lstsq 将返回该精确解。

步骤：

豆包AI编程

豆包推出的AI编程助手

483

查看详情

将约束表示为矩阵形式 C X = d。 对于给定的约束：
- 0.5 * y1 + 0.5 * y2 = 0
- 0.5 * x3 + 0.5 * x4 = 0
- 0.5 * y3 + 0.5 * y4 = 0
我们可以构建一个 3x8 的矩阵 C 和一个 3x1 的向量 d。 X = [x1, y1, x2, y2, x3, y3, x4, y4]

C 矩阵的行对应每个约束，列对应 X 中的变量： C = [[0, 0.5, 0, 0.5, 0, 0, 0, 0],[0, 0, 0, 0, 0.5, 0, 0.5, 0],[0, 0, 0, 0, 0, 0.5, 0, 0.5]]d = [[0], [0], [0]]
增广 A 和 b。 将 C 垂直堆叠到 A 下方，将 d 垂直堆叠到 b 下方。

A_aug = np.vstack([A, C])b_aug = np.vstack([b, d])
使用 np.linalg.lstsq 求解。np.linalg.lstsq(A_aug, b_aug, rcond=None) 将返回增广系统的最小二乘解。

# 原始 A 和 b (与上文相同)
# A = ...
# b = ...

# 1. 构建约束矩阵 AC 和约束向量 bC
AC = np.zeros([3, A.shape[1]]) # 3个约束，8个变量
bC = np.zeros((3, 1))

# 填充 AC 矩阵
# X = [x1, y1, x2, y2, x3, y3, x4, y4]
# 索引: x[0]=x1, x[1]=y1, x[2]=x2, x[3]=y2, x[4]=x3, x[5]=y3, x[6]=x4, x[7]=y4

# 约束 1: 0.5*(y1 + y2) = 0  => 0.5*x[1] + 0.5*x[3] = 0
AC[0][[1, 3]] = 0.5

# 约束 2: 0.5*(x3 + x4) = 0  => 0.5*x[4] + 0.5*x[6] = 0
AC[1][[4, 6]] = 0.5

# 约束 3: 0.5*(y3 + y4) = 0  => 0.5*x[5] + 0.5*x[7] = 0
AC[2][[5, 7]] = 0.5

# bC 向量已初始化为零

# 2. 增广系统
A_augmented = np.vstack([A, AC])
b_augmented = np.vstack([b, bC])

print("增广后的 A 矩阵形状:", A_augmented.shape)
print("增广后的 b 向量形状:", b_augmented.shape)

# 3. 使用 np.linalg.lstsq 求解增广系统
# rcond=None 禁用 rcond 警告
x_lstsq, residuals, rank, singular_values = np.linalg.lstsq(A_augmented, b_augmented, rcond=None)

print("\nnp.linalg.lstsq 找到的解 X:")
print(x_lstsq.flatten())

# 验证约束条件
print("\n验证约束条件 (应接近于0):")
# 注意：x_lstsq 是一个列向量，需要展平或适当索引
print(np.dot(AC, x_lstsq).flatten())

# 验证原始 AX 与 b 的匹配程度
print("\n验证原始 AX 与 b 的匹配程度:")
print(np.matmul(A, x_lstsq).flatten())
print("\n期望的 b 向量 (原始):")
print(b.flatten())

# 检查原始 AX 和 b 之间的残差
original_residuals = np.matmul(A, x_lstsq) - b
print("\n原始 AX 与 b 的残差:")
print(original_residuals.flatten())
print("原始 AX 与 b 的残差平方和:", np.sum(original_residuals**2))

登录后复制

通过这种方法，np.linalg.lstsq 会找到一个 X，它在最小二乘意义上最佳地满足了所有 11 个方程（8个原始方程 + 3个约束方程）。这意味着它会尽可能地使 AX 接近 b，同时精确地满足（或尽可能接近满足，如果系统高度不一致）所有线性约束。

在上述示例中，np.linalg.lstsq 找到的解 x_lstsq 在满足约束的同时，会使 np.matmul(A, x_lstsq) 的结果更接近原始 b 向量，或者在整个增广系统上达到最小的残差。

总结与注意事项

选择方法：
- 当约束是线性的，并且目标是找到一个既满足原始方程组（或其最佳最小二乘近似）又满足所有约束的解时，增广系统结合 np.linalg.lstsq 是首选。这种方法直接将约束融入到待解系统中，求得的解会同时考虑所有条件。
- 当约束是非线性的，或者目标函数除了 ||AX - b||^2 之外还有其他复杂的项需要最小化时，scipy.optimize.minimize 及其变体（如 SLSQP）是更通用的选择。但需注意，minimize 仅保证目标函数在约束下的最小值，不一定能使 AX 精确等于 b。
过定系统： 当增广后的方程数多于未知数时（如本例中的 11 个方程，8 个未知数），系统是过定的。np.linalg.lstsq 能够稳健地处理过定系统，找到一个最小二乘意义上的最佳拟合解。
rcond 参数： np.linalg.lstsq 中的 rcond 参数用于控制小奇异值的处理，以防止在病态矩阵情况下产生不稳定的解。在较新版本的 NumPy 中，推荐将其设置为 None 以使用默认行为。