在Python中为16个变量拟合17个不等式约束的求解问题

我有16个变量（x[0]到x[15]），需要找到满足以下17个不等式约束的取值：

约束条件 #

1:  (x[0] + x[1] + x[2] + x[4] + x[6] + x[7])/(x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + x[7]) >= 0.4
2:  (x[3] + x[5])/(x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + x[7]) <= 0.6
3:  x[0]/(x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + x[7]) <= 0.1
4:  (x[0] + x[2] + x[4] + x[6])/(x[0] + x[1] + x[2] + x[3] + x[4] + x[5] + x[6] + x[7]) <= 0.1
5:  (1.2*x[11]*x[3] + 1.2*x[13]*x[5])/(x[3] + x[5]) <= 520
6:  (1.2*x[11]*x[3] + 1.2*x[13]*x[5])/(x[3] + x[5]) >= 470
7:  (x[0]*x[8] + x[1]*x[9] + x[10]*x[2] + x[12]*x[4] + x[14]*x[6] + x[15]*x[7])/(x[0] + x[1] + x[2] + x[4] + x[6] + x[7]) <= 420
8:  (x[3] + x[7])/x[5] >= 0.05
9:  (x[3] + x[7])/x[5] <= 0.2
10:  x[2]/(x[3] + x[7]) >= 0.05
11:  x[2]/(x[3] + x[7]) <= 0.15
12:  x[5]/(x[3] + x[5]) >= 0.95
13:  0.833333333333333/x[11] >= 376
14:  0.833333333333333/x[11] <= 424 
15:  x[13]/x[11] >= 0.7
16:  x[13]/x[11] <= 0.82
17:  1.2*x[11]*x[3] + 1.2*x[13]*x[5] <= 317300.0

变量范围 #

• x[0]到x[7]：20 ≤ 取值 ≤ 600
• x[8]到x[15]：200 ≤ 取值 ≤ 600

过往尝试 #

曾用scipy.optimize.minimize()的SLSQP方法，以变量之和为目标函数求解，但部分约束未被满足。无需最小化变量值，仅需找到符合所有约束的可行解。

1. 约束预处理：消除分式 + 简化逻辑 #

所有变量下限为正数，分母恒大于0，可直接将分式约束转化为非分式形式，避免求解器因分母波动出错，同时简化约束逻辑：

• 约束1 → x[0]+x[1]+x[2]+x[4]+x[6]+x[7] ≥ 0.4×sum(x[0:8])
• 约束2 → x[3]+x[5] ≤ 0.6×sum(x[0:8])
• 约束3 → x[0] ≤ 0.1×sum(x[0:8])
• 约束4 → x[0]+x[2]+x[4]+x[6] ≤ 0.1×sum(x[0:8])
• 约束5 → 1.2x[11]x[3] + 1.2x[13]x[5] ≤ 520(x[3]+x[5])
• 约束6 → 1.2x[11]x[3] + 1.2x[13]x[5] ≥ 470(x[3]+x[5])
• 约束7 → x[0]x[8]+x[1]x[9]+x[10]x[2]+x[12]x[4]+x[14]x[6]+x[15]x[7] ≤ 420(x[0]+x[1]+x[2]+x[4]+x[6]+x[7])
• 约束8 → x[3]+x[7] ≥ 0.05x[5]
• 约束9 → x[3]+x[7] ≤ 0.2x[5]
• 约束10 → x[2] ≥ 0.05(x[3]+x[7])
• 约束11 → x[2] ≤ 0.15(x[3]+x[7])
• 约束12 → 简化为x[5] ≥ 19x[3]（原约束变形后推导得出）
• 约束13 → 0.833333333333333 ≥ 376x[11]
• 约束14 → 424x[11] ≥ 0.833333333333333
• 约束15 → x[13] ≥ 0.7x[11]
• 约束16 → x[13] ≤ 0.82x[11]
• 约束17 → 保持原形式

2. 调整scipy求解策略：聚焦可行解 #

SLSQP对初始点敏感度高，调整目标函数与初始点，让求解器专注于寻找可行解：

import numpy as np
from scipy.optimize import minimize

# 生成符合部分约束的初始点：满足x5≥19x3，变量取范围内合理值
x0 = np.zeros(16)
x0[0:3] = [300, 300, 300]
x0[3] = 25  # 满足x3≤600/19≈31.57
x0[5] = 19*25  # 475，符合20-600范围
x0[4] = 300
x0[6:8] = [300, 300]
x0[8:16] = [400]*8  # x8-x15取中间值

# 目标函数设为常数0，仅寻找可行解
def objective(x):
    return 0.0

# 定义转化后的约束
constraints = [
    {'type': 'ineq', 'fun': lambda x: (x[0]+x[1]+x[2]+x[4]+x[6]+x[7]) - 0.4*sum(x[0:8])},
    {'type': 'ineq', 'fun': lambda x: 0.6*sum(x[0:8]) - (x[3]+x[5])},
    {'type': 'ineq', 'fun': lambda x: 0.1*sum(x[0:8]) - x[0]},
    {'type': 'ineq', 'fun': lambda x: 0.1*sum(x[0:8]) - (x[0]+x[2]+x[4]+x[6])},
    {'type': 'ineq', 'fun': lambda x: 520*(x[3]+x[5]) - (1.2*x[11]*x[3] + 1.2*x[13]*x[5])},
    {'type': 'ineq', 'fun': lambda x: (1.2*x[11]*x[3] + 1.2*x[13]*x[5]) - 470*(x[3]+x[5])},
    {'type': 'ineq', 'fun': lambda x: 420*(x[0]+x[1]+x[2]+x[4]+x[6]+x[7]) - (x[0]*x[8]+x[1]*x[9]+x[10]*x[2]+x[12]*x[4]+x[14]*x[6]+x[15]*x[7])},
    {'type': 'ineq', 'fun': lambda x: (x[3]+x[7]) - 0.05*x[5]},
    {'type': 'ineq', 'fun': lambda x: 0.2*x[5] - (x[3]+x[7])},
    {'type': 'ineq', 'fun': lambda x: x[2] - 0.05*(x[3]+x[7])},
    {'type': 'ineq', 'fun': lambda x: 0.15*(x[3]+x[7]) - x[2]},
    {'type': 'ineq', 'fun': lambda x: x[5] - 19*x[3]},
    {'type': 'ineq', 'fun': lambda x: 0.833333333333333 - 376*x[11]},
    {'type': 'ineq', 'fun': lambda x: 424*x[11] - 0.833333333333333},
    {'type': 'ineq', 'fun': lambda x: x[13] - 0.7*x[11]},
    {'type': 'ineq', 'fun': lambda x: 0.82*x[11] - x[13]},
    {'type': 'ineq', 'fun': lambda x: 317300.0 - (1.2*x[11]*x[3] + 1.2*x[13]*x[5])}
]

# 设置变量边界
bounds = [(20,600)]*8 + [(200,600)]*8

# 求解：增加迭代次数与精度
result = minimize(objective, x0, method='SLSQP', bounds=bounds, constraints=constraints,
                  options={'maxiter': 1000, 'ftol': 1e-8, 'disp': True})

# 输出结果与约束验证
print("求解成功:", result.success)
if result.success:
    print("\n可行解：")
    for idx, val in enumerate(result.x):
        print(f"x[{idx}] = {val:.4f}")
    print("\n约束验证（值≥-1e-6即满足）：")
    for i, con in enumerate(constraints):
        val = con['fun'](result.x)
        print(f"约束{i+1}：{val:.4f} → {'满足' if val >= -1e-6 else '不满足'}")

3. 备选方案：用Pyomo建模求解 #

如果scipy仍无法收敛，可使用Pyomo这类专业优化库，搭配IPOPT等求解器，更适合复杂约束的可行性问题：

from pyomo.environ import *

# 创建模型
model = ConcreteModel()

# 定义变量与边界
model.x = Var(range(16), domain=Reals)
for i in range(8):
    model.x[i].setlb(20)
    model.x[i].setub(600)
for i in range(8,16):
    model.x[i].setlb(200)
    model.x[i].setub(600)

# 添加约束
model.con1 = Constraint(expr=(model.x[0]+model.x[1]+model.x[2]+model.x[4]+model.x[6]+model.x[7]) >= 0.4*sum(model.x[i] for i in range(8)))
model.con2 = Constraint(expr=(model.x[3]+model.x[5]) <= 0.6*sum(model.x[i] for i in range(8)))
model.con3 = Constraint(expr=model.x[0] <= 0.1*sum(model.x[i] for i in range(8)))
model.con4 = Constraint(expr=(model.x[0]+model.x[2]+model.x[4]+model.x[6]) <= 0.1*sum(model.x[i] for i in range(8)))
model.con5 = Constraint(expr=(1.2*model.x[11]*model.x[3] + 1.2*model.x[13]*model.x[5]) <= 520*(model.x[3]+model.x[5]))
model.con6 = Constraint(expr=(1.2*model.x[11]*model.x[3] + 1.2*model.x[13]*model.x[5]) >= 470*(model.x[3]+model.x[5]))
model.con7 = Constraint(expr=(model.x[0]*model.x[8]+model.x[1]*model.x[9]+model.x[10]*model.x[2]+model.x[12]*model.x[4]+model.x[14]*model.x[6]+model.x[15]*model.x[7]) <= 420*(model.x[0]+model.x[1]+model.x[2]+model.x[4]+model.x[6]+model.x[7]))
model.con8 = Constraint(expr=(model.x[3]+model.x[7]) >= 0.05*model.x[5])
model.con9 = Constraint(expr=(model.x[3]+model.x[7]) <= 0.2*model.x[5])
model.con10 = Constraint(expr=model.x[2] >= 0.05*(model.x[3]+model.x[7]))
model.con11 = Constraint(expr=model.x[2] <= 0.15*(model.x[3]+model.x[7]))
model.con12 = Constraint(expr=model.x[5] >= 19*model.x[3])
model.con13 = Constraint(expr=0.833333333333333 >= 376*model.x[11])
model.con14 = Constraint(expr=424*model.x[11] >= 0.833333333333333)
model.con15 = Constraint(expr=model.x[13] >= 0.7*model.x[11])
model.con16 = Constraint(expr=model.x[13] <= 0.82*model.x[11])
model.con17 = Constraint(expr=1.2*model.x[11]*model.x[3] + 1.2*model.x[13]*model.x[5] <= 317300.0)

# 目标函数设为常数0
model.obj = Objective(expr=0.0, sense=minimize)

# 求解（需提前安装IPOPT求解器）
solver = SolverFactory('ipopt')
result = solver.solve(model, tee=True)

# 输出结果
print("\n求解状态:", result.solver.status)
if result.solver.status == 'ok':
    print("\n可行解：")
    for i in range(16):
        print(f"x[{i}] = {model.x[i].value:.4f}")

关键提示 #

• 约束12限定x[3]最大约为31.57（因x[5]≤600），可在变量边界中直接设置x[3]上限为32，缩小搜索范围，提升求解效率。
• 若多次求解仍失败，需检查约束是否存在矛盾（比如多个约束是否无法同时满足）。

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