scipy.optimize.linprog未返回预期最优解问题咨询

问题描述 #

使用Scipy的linprog求解线性规划问题时，已知该函数默认求最小值，将目标函数乘以-1可转化为求最大值，且官方示例验证有效。但在自身问题中：

• 最小化目标函数[0.5, 0.5]时，得到解[0.61538462 0.38461538]，符合预期；
• 最大化目标函数（即最小化[-0.5, -0.5]）时，预期得到解[0.0, 1.0]，但实际结果仍为[0.61538462 0.38461538]。

测试代码：

from scipy.optimize import linprog

obj_fct1 = [0.5, 0.5]
obj_fct2 = [-0.5, -0.5]
lhs_ineq = [[1.5, 0.2]]
rhs_ineq = [1]
lhs_eq = [[1,1]]
rhs_eq = [1]
bnds = [(0, 1),
          (0, 1)]

opt_sol1 = linprog(c=obj_fct1,
                  A_ub=lhs_ineq,
                  b_ub=rhs_ineq,
                  A_eq=lhs_eq,
                  b_eq=rhs_eq,
                  bounds=bnds)
print(opt_sol1.x)
print("------------------------------------------")
opt_sol2 = linprog(c=obj_fct2,
                  A_ub=lhs_ineq,
                  b_ub=rhs_ineq,
                  A_eq=lhs_eq,
                  b_eq=rhs_eq,
                  bounds=bnds)
print(opt_sol2.x)

# 输出：
# [0.61538462 0.38461538]
# ------------------------------------------
# [0.61538462 0.38461538]

核心原因：目标函数在约束下为常数 #

你的等式约束是x1 + x2 = 1，代入目标函数0.5x1 + 0.5x2可得：
0.5x1 + 0.5x2 = 0.5*(x1+x2) = 0.5*1 = 0.5
也就是说，所有满足约束条件的可行解，对应的目标函数值都是0.5。因此：

• 最小化0.5x1+0.5x2时，所有可行解都是最优解；
• 最大化0.5x1+0.5x2（即最小化-0.5x1-0.5x2）时，目标函数值恒为-0.5，所有可行解也都是最优解。

linprog在存在多个最优解时，只会返回求解过程中找到的第一个可行最优解，而非所有最优解，因此两次求解返回了同一个结果。

另外验证[0.0,1.0]确实是可行解：代入不等式约束1.5*0 + 0.2*1 = 0.2 ≤ 1，满足；等式约束0+1=1，满足；边界0≤x1≤1, 0≤x2≤1，满足。所以它是最优解之一。

如何获取预期的[0.0,1.0]解 #

因为所有可行解都是最优解，你可以通过以下方式引导求解器返回目标点：

1. 添加微小扰动到目标函数 #

微调最大化的目标函数，让求解器偏向x2更大的方向：

# 微调目标函数，让x2的权重略大于x1
obj_fct2_perturbed = [-0.5, -0.5001]
opt_sol2 = linprog(c=obj_fct2_perturbed,
                   A_ub=lhs_ineq,
                   b_ub=rhs_ineq,
                   A_eq=lhs_eq,
                   b_eq=rhs_eq,
                   bounds=bnds)
print(opt_sol2.x)  # 输出接近[0.0, 1.0]

2. 指定初始点（适用于支持该参数的求解器） #

使用method='highs'求解器，指定初始点为[0.0,1.0]：

opt_sol2 = linprog(c=obj_fct2,
                   A_ub=lhs_ineq,
                   b_ub=rhs_ineq,
                   A_eq=lhs_eq,
                   b_eq=rhs_eq,
                   bounds=bnds,
                   method='highs',
                   x0=[0.0, 1.0])
print(opt_sol2.x)  # 输出[0.0, 1.0]

3. 手动取用边界可行解 #

由于问题变量少，可直接计算所有边界可行解：

• 当x1=0时，x2=1，满足所有约束；
• 当x1≈0.61538时，x2≈0.38462，是不等式约束的边界点。
  这些点都是最优解，可直接取用。

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