如何在Scipy中解决等式约束超独立变量的最小化问题（含倍数约束）

报错原因 #

你添加的约束eq_cons = {'type': 'eq', 'fun': lambda x: x%1}会返回一个和输入x同维度的数组（4个元素），这相当于一次性添加了4个等式约束（要求每个x_i的小数部分为0）。加上原有的3个线性等式约束，总共7个等式约束，但优化变量只有4个，远超变量自由度，因此SLSQP报错"More equality constraints than independent variables"。

此外，x%1这类离散模运算约束本质是非光滑、非连续的，SLSQP作为连续优化器，本身就不适合处理这类整数/倍数约束。

正确解决思路 #

要实现"x为常数k的倍数"的约束，需要使用**混合整数线性规划（MILP）**工具，因为你的目标函数和约束都可以转化为线性形式，适合用scipy.optimize.milp求解（Scipy 1.9+版本支持）。

步骤1：问题转化 #

假设要求每个x_i是常数k的倍数，令x = k * y，其中y_i为整数。这样问题就转化为对整数变量y的线性规划：

• 线性约束：A*(k*y) = B → (k*A)y = B
• 变量边界：ceil(l_i/k) ≤ y_i ≤ floor(u_i/k)（确保x_i在原边界内）
• 目标函数：abs(sum(x - (u+l)/2))可以转化为线性形式：引入辅助变量s，最小化s，同时满足s ≥ sum(x - mid)和s ≥ sum(mid - x)（其中mid=(u+l)/2）

代码示例（以k=1为例，即x为整数） #

import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds

# 原问题参数
A = np.array([[0.106667, 0.1333, 0.1333, 0.01],
              [0.02, 0.6667, 0.1333, 0.12],
              [0.0933, 0.06667, 0.6, 0.01]])
B = np.array([27, 57, 28])  # 转为一维数组简化计算
l = np.array([100, 40, 10, 50])
u = np.array([200, 80, 20, 150])
mid = (u + l) / 2
sum_mid = np.sum(mid)
k = 1  # 自定义倍数常数

# 变量定义：y(4维整数) + s(1维连续辅助变量)，共5个变量
n_vars = 4 + 1
# 目标函数：最小化s，系数数组前4位对应y，最后一位对应s
c = np.zeros(n_vars)
c[-1] = 1

# 约束1：原线性等式约束 A*x = B → A*(k*y) = B
A_eq = k * A
lin_con1 = LinearConstraint(A_eq, B, B)

# 约束2：sum(x) - sum_mid ≤ s → sum(k*y) - s ≤ sum_mid
coeffs2 = np.hstack([k*np.ones(4), -1])
lin_con2 = LinearConstraint(coeffs2, -np.inf, sum_mid)

# 约束3：sum_mid - sum(x) ≤ s → -sum(k*y) - s ≤ -sum_mid
coeffs3 = np.hstack([-k*np.ones(4), -1])
lin_con3 = LinearConstraint(coeffs3, -np.inf, -sum_mid)

# 变量边界：y的边界由原x边界转换而来，s≥0
y_lower = np.ceil(l / k)
y_upper = np.floor(u / k)
bounds = Bounds(np.hstack([y_lower, 0]), np.hstack([y_upper, np.inf]))

# 指定整数变量：前4个变量(y)是整数，最后一个(s)是连续变量
integrality = np.array([1, 1, 1, 1, 0])

# 求解
res = milp(c=c, constraints=[lin_con1, lin_con2, lin_con3], bounds=bounds, integrality=integrality)

# 输出结果
print("优化状态:", res.status)
print("最优目标值:", res.fun)
# 还原x值：x = k*y
x_opt = k * res.x[:4]
print("最优x值:", x_opt)
# 验证线性约束是否满足
print("线性约束误差:", A.dot(x_opt) - B)

适配任意k值 #

如果需要x是其他常数k的倍数（比如k=5），只需修改代码中的k值即可，其余逻辑自动适配：

• A_eq会自动变为5*A
• y的边界会自动计算为ceil(l/5)和floor(u/5)
• 目标函数的约束系数也会自动调整为5*np.ones(4)

备选方案（启发式） #

如果不想使用MILP，可以先求解无整数约束的连续最优解，再将x调整为最近的k的倍数，最后验证是否满足线性约束。但这种方法无法保证找到可行解，仅适合对精度要求不高的场景：

# 先求解连续解
from scipy.optimize import minimize, Bounds

def objfun(x):
    return abs(np.sum(x - (u+l)/2))

x0 = np.array([150, 60, 15, 100])
eq_cons1 = {'type': 'eq', 'fun': lambda x: np.matmul(A[0,:],x)-B[0]}
eq_cons2 = {'type': 'eq', 'fun': lambda x: np.matmul(A[1,:],x)-B[1]}
eq_cons3 = {'type': 'eq', 'fun': lambda x: np.matmul(A[2,:],x)-B[2]}

res_cont = minimize(objfun, x0, method='SLSQP',
                    constraints=[eq_cons1, eq_cons2, eq_cons3],
                    bounds=Bounds(l, u))

# 调整为k的倍数
k = 1
x_adj = np.round(res_cont.x / k) * k
# 裁剪到边界内
x_adj = np.clip(x_adj, l, u)

# 验证约束
print("调整后x值:", x_adj)
print("线性约束误差:", A.dot(x_adj) - B)

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