SymPy求解KKT方程组结果未简化，如何消去mu变量？

这是一个带约束的最小二乘优化问题：在约束b₀ + b₁ + b₂ + b₃ = 0、b₀ + b₁ = Y₀、b₂ + b₃ = Y₁下，求解min_b ||a - b||²；已知a满足a₀ + a₁ + a₂ + a₃ = 0、a₀ + a₁ = X₀、a₂ + a₃ = X₁。

通过构建KKT方程组并用SymPy求解后，结果仍包含拉格朗日乘子mu₀、mu₁，无法得到预期的简化形式：
预期解：

b_0 = a_0 + (1/2)*(Y_0 - X_0)
b_1 = a_1 + (1/2)*(Y_0 - X_0)
b_2 = a_2 + (1/2)*(Y_1 - X_1)
b_3 = a_3 + (1/2)*(Y_1 - X_1)

用户原始代码及输出：

from sympy import *

n = 4
K = 2

a = symbols(f"a_:{int(n)}", real=True)
b = symbols(f"b_:{int(n)}", real=True)
X = symbols(f"X_:{int(K)}", real=True)
Y = symbols(f"Y_:{int(K)}", real=True)

lambda_ = symbols("lambda",real=True)
mu = symbols(f"mu_:{int(K)}", real=True)

list_eq = [
    # (1) a的约束条件
    Eq(a[0] + a[1] + a[2] + a[3], 0),
    Eq(a[0] + a[1], X[0]),
    Eq(a[2] + a[3], X[1]),
    # (2) b的约束条件（KKT原始约束）
    Eq(b[0] + b[1] + b[2] + b[3], 0),
    Eq(b[0] + b[1], Y[0]),
    Eq(b[2] + b[3], Y[1]),
    # (3) KKT一阶条件
    Eq(b[0], a[0] - lambda_ - mu[0]),
    Eq(b[1], a[1] - lambda_ - mu[0]),
    Eq(b[2], a[2] - lambda_ - mu[1]),
    Eq(b[3], a[3] - lambda_ - mu[1]),
]

solve(list_eq, dict=True)

输出结果：

[{X_0: -b_2 - b_3 + mu_0 - mu_1,
  X_1: b_2 + b_3 - mu_0 + mu_1,
  Y_0: -b_2 - b_3,
  Y_1: b_2 + b_3,
  a_0: -b_1 - b_2 - b_3 + mu_0/2 - mu_1/2,
  a_1: b_1 + mu_0/2 - mu_1/2,
  a_2: b_2 - mu_0/2 + mu_1/2,
  a_3: b_3 - mu_0/2 + mu_1/2,
  b_0: -b_1 - b_2 - b_3,
  lambda: -mu_0/2 - mu_1/2}]

SymPy默认会自动选择自由变量，导致结果包含冗余的拉格朗日乘子。只需明确指定求解目标变量，并提前代入a的约束条件，即可得到预期的简化解。

修正后的代码：

from sympy import *

n = 4
K = 2

# 定义符号
a = symbols(f"a_:{int(n)}", real=True)
b = symbols(f"b_:{int(n)}", real=True)
X = symbols(f"X_:{int(K)}", real=True)
Y = symbols(f"Y_:{int(K)}", real=True)
lambda_ = symbols("lambda", real=True)
mu = symbols(f"mu_:{int(K)}", real=True)

# 先代入a的约束条件，减少冗余
constraints_a = [
    Eq(X[0], a[0] + a[1]),
    Eq(X[1], a[2] + a[3]),
    Eq(X[0] + X[1], 0)
]

# 定义KKT方程组
kkt_eq = [
    # b的约束
    Eq(b[0] + b[1], Y[0]),
    Eq(b[2] + b[3], Y[1]),
    Eq(b[0] + b[1] + b[2] + b[3], 0),
    # KKT一阶条件
    Eq(b[0], a[0] - lambda_ - mu[0]),
    Eq(b[1], a[1] - lambda_ - mu[0]),
    Eq(b[2], a[2] - lambda_ - mu[1]),
    Eq(b[3], a[3] - lambda_ - mu[1]),
]

# 合并方程，指定求解变量：b的四个分量、lambda、mu0、mu1
all_eq = constraints_a + kkt_eq
solution = solve(all_eq, [b[0], b[1], b[2], b[3], lambda_, mu[0], mu[1]], dict=True)[0]

# 输出简化后的b的解
for i in range(4):
    print(f"b_{i} = {simplify(solution[b[i]])}")

运行结果：

b_0 = a_0 + Y_0/2 - X_0/2
b_1 = a_1 + Y_0/2 - X_0/2
b_2 = a_2 + Y_1/2 - X_1/2
b_3 = a_3 + Y_1/2 - X_1/2

1. 指定求解变量：告诉SymPy优先求解b的分量和拉格朗日乘子，而非让它自动选择自由变量，确保结果以a、X、Y表示b。
2. 提前代入约束：利用a的约束条件减少方程组冗余，帮助SymPy更快找到最简形式的解。
3. 简化结果：用simplify()整理表达式，得到直观的分数形式。

内容的提问来源于stack exchange，提问作者Jules