Python线性/非线性规划优化遇阻：Scipy/Pulp求解失败，求可行方法

首先澄清一下约束表述的歧义：从你的代码逻辑来看，你应该是想表达至少有2个$j$满足$\frac{x_j a_j}{\sum x_i a_i} \geq 0.15$、至少有3个$j$满足$\frac{x_j a_j}{\sum x_i a_i} \geq 0.07$——如果是“任意两个/三个$j$”的话，会出现约束矛盾（比如当$n=9$时，9个比例都≥0.15的话总和会超过1），所以下面的方案基于这个合理的理解。

问题回顾 #

你需要求解的优化问题：

• 目标函数：最大化 $\sum_{i=1}^n x_i a_i$（其中$a_i$是给定的share数组，$1 \leq x_i \leq 100$且为整数）
• 约束条件：
  1. 所有$j$的比例 $\frac{x_j a_j}{\sum x_i a_i} \leq 0.35$
  2. 至少2个$j$的比例 $\geq 0.15$
  3. 至少3个$j$的比例 $\geq 0.07$

现有Scipy代码的问题 #

你的SLSQP求解失败主要有几个核心原因：

1. 非光滑约束函数：你用计数统计满足/不满足约束的数量来构建约束，这种离散的计数函数是不连续、不可微的，而SLSQP这类基于梯度的优化器只能处理连续可微的约束。
2. 约束逻辑错误：比如constraint2里多乘了100，导致约束条件和你的原始需求不符。
3. 初始点问题：全100的初始点让所有比例完全相同，此时优化器的梯度矩阵出现奇异，无法继续迭代。

正确方案：转化为混合整数线性规划（MILP） #

这个问题可以通过变量替换转化为线性约束的形式，再引入少量0-1整数变量来处理“至少k个满足”的约束，这样就能用成熟的MILP求解器（比如PuLP搭配CBC）来稳定求解。

变量替换思路 #

令 $S = \sum_{i=1}^n x_i a_i$（也就是我们要最大化的目标），再令 $w_i = x_i a_i$，那么：

• $w_i$ 的范围是 $a_i \leq w_i \leq 100a_i$（因为$1 \leq x_i \leq 100$）
• 目标变为最大化 $S$
• 所有比例约束可以转化为关于$w_i$和$S$的线性约束

对于“至少k个满足”的约束，我们引入0-1变量$y_j$（$y_j=1$表示第j个比例满足要求，$y_j=0$则不满足），通过线性约束强制当$y_j=1$时，$w_j$必须符合比例要求，同时要求$\sum y_j \geq k$。

PuLP实现代码 #

import pulp

# 给定的a_i数组（对应你的share）
share = [1595798.061003, 1595798.061003, 1595798.061003, 1595798.061003,
         6335021.83000001, 6335021.83000001, 6335021.83000001, 6335021.83000001,
         42842994.4958]
n = len(share)

# 创建最大化问题实例
prob = pulp.LpProblem("Maximize_Weighted_Sum", pulp.LpMaximize)

# 定义变量：
# x_i 是整数变量，范围[1,100]
x = [pulp.LpVariable(f"x_{i}", lowBound=1, upBound=100, cat='Integer') for i in range(n)]
# S是目标函数值，即sum(x_i * a_i)
S = pulp.LpVariable("S", lowBound=0)
# y_j是0-1变量，标记是否满足比例≥0.15
y = [pulp.LpVariable(f"y_{i}", cat='Binary') for i in range(n)]
# z_j是0-1变量，标记是否满足比例≥0.07
z = [pulp.LpVariable(f"z_{i}", cat='Binary') for i in range(n)]

# 目标函数：最大化S
prob += S

# 约束：S等于所有x_i*a_i的和
prob += pulp.lpSum([x[i] * share[i] for i in range(n)]) == S

# 约束1：每个比例≤0.35 → x_j*a_j ≤ 0.35*S
for j in range(n):
    prob += x[j] * share[j] <= 0.35 * S

# 约束2：至少2个比例≥0.15
prob += pulp.lpSum(y) >= 2
for j in range(n):
    # 当y_j=1时，强制x_j*a_j ≥ 0.15*S；y_j=0时，约束自动满足
    prob += x[j] * share[j] >= 0.15 * S * y[j]

# 约束3：至少3个比例≥0.07
prob += pulp.lpSum(z) >= 3
for j in range(n):
    prob += x[j] * share[j] >= 0.07 * S * z[j]

# 求解（关闭日志输出，若需要调试可去掉msg=0）
prob.solve(pulp.PULP_CBC_CMD(msg=0))

# 输出结果
print("求解状态:", pulp.LpStatus[prob.status])
print("最大化的目标值S:", pulp.value(S))
print("\n各变量结果:")
for i in range(n):
    ratio = (pulp.value(x[i]) * share[i]) / pulp.value(S)
    print(f"第{i+1}个变量: x_i={pulp.value(x[i])}, 比例={ratio:.4f}")

为什么这个方法有效？ #

PuLP使用的CBC求解器是专门为混合整数线性规划设计的，能够处理我们引入的0-1变量，同时所有约束都转化为了线性形式，避免了原代码中的非线性和非光滑问题，求解稳定性和可靠性都远高于用SLSQP处理非光滑约束的方式。

其他可选方案 #

如果你不想用整数变量，也可以考虑使用支持非线性规划的求解器（比如IPOPT，可通过Pyomo或CasADi调用），但需要将约束改写为连续可微的形式，不过这种方法的稳定性不如MILP。

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