投资组合优化MILP问题：如何在CVXPY中实现剩余项约束？

问题背景 #

这是一个投资组合优化领域的混合整数线性规划（MILP）问题，使用CVXPY工具实现。现有代码已完成约束1和约束2的实现，需要添加约束3：

• 约束3定义：others为[b, c, d]中除去b以及其中最大的两个元素对应的sum值后的剩余元素。若others非空，其中每个元素i需满足sum(x*i) ≤10。
  • 示例1：当[sum(x*b), sum(x*c), sum(x*d)] = [1,2,3]时，最大两个sum对应c、d，others为空，无需约束；
  • 示例2：当数组为[3,2,1]时，最大两个sum对应b、c，others为d，需满足sum(x*d) ≤10。

原代码（修正后） #

import numpy as np
import cvxpy as cp

n = 10

a = np.random.randint(1, 10, size=n)
b = np.random.randint(1, 10, size=n)
c = np.random.randint(1, 10, size=n)
d = np.random.randint(1, 10, size=n)

x = cp.Variable(shape=n, boolean=True)

# 目标函数
objective = cp.Maximize(cp.sum(cp.multiply(x, a)))

# 约束
constraints = []
# 约束1：修正原代码括号错误
constraints.append(cp.sum(cp.multiply(x, b)) <= 50)
# 约束2：最大两个sum的和 ≤100
S_b = cp.sum(cp.multiply(x, b))
S_c = cp.sum(cp.multiply(x, c))
S_d = cp.sum(cp.multiply(x, d))
constraints.append(cp.sum_largest(cp.hstack([S_b, S_c, S_d]), 2) <= 100)

约束3的实现 #

由于CVXPY无法直接处理条件判断逻辑，我们通过二进制辅助变量+大M法将约束3转化为线性约束：

步骤1：定义辅助变量与大M值 #

# 定义大M：取sum的最大可能值（每个元素最大9，n=10，故最大sum=90，M取100足够）
M = 100

# 二进制辅助变量：y_c=1表示S_c是三个sum中的最小值；y_d=1表示S_d是三个sum中的最小值
y_c = cp.Variable(boolean=True)
y_d = cp.Variable(boolean=True)

# 三个sum中最多有一个最小值（若最小值是b，则y_c和y_d都为0）
constraints.append(y_c + y_d <= 1)

步骤2：用大M约束定义辅助变量的逻辑 #

# 当y_c=1时，S_c ≤ S_b且S_c ≤ S_d（即S_c是最小值）
constraints.append(S_c <= S_b + M * (1 - y_c))
constraints.append(S_c <= S_d + M * (1 - y_c))

# 当y_d=1时，S_d ≤ S_b且S_d ≤ S_c（即S_d是最小值）
constraints.append(S_d <= S_b + M * (1 - y_d))
constraints.append(S_d <= S_c + M * (1 - y_d))

步骤3：添加约束3的核心要求 #

# 若S_c是最小值（且不是b），则S_c ≤10；否则约束自动放宽
constraints.append(S_c <= 10 + M * (1 - y_c))
# 若S_d是最小值（且不是b），则S_d ≤10；否则约束自动放宽
constraints.append(S_d <= 10 + M * (1 - y_d))

完整求解代码 #

import numpy as np
import cvxpy as cp

n = 10

a = np.random.randint(1, 10, size=n)
b = np.random.randint(1, 10, size=n)
c = np.random.randint(1, 10, size=n)
d = np.random.randint(1, 10, size=n)

x = cp.Variable(shape=n, boolean=True)

# 目标函数
objective = cp.Maximize(cp.sum(cp.multiply(x, a)))

# 约束
constraints = []
# 约束1
S_b = cp.sum(cp.multiply(x, b))
constraints.append(S_b <= 50)
# 约束2
S_c = cp.sum(cp.multiply(x, c))
S_d = cp.sum(cp.multiply(x, d))
constraints.append(cp.sum_largest(cp.hstack([S_b, S_c, S_d]), 2) <= 100)

# 约束3实现
M = 100
y_c = cp.Variable(boolean=True)
y_d = cp.Variable(boolean=True)
constraints.append(y_c + y_d <= 1)

# 定义辅助变量逻辑
constraints.append(S_c <= S_b + M * (1 - y_c))
constraints.append(S_c <= S_d + M * (1 - y_c))
constraints.append(S_d <= S_b + M * (1 - y_d))
constraints.append(S_d <= S_c + M * (1 - y_d))

# 核心约束
constraints.append(S_c <= 10 + M * (1 - y_c))
constraints.append(S_d <= 10 + M * (1 - y_d))

# 求解模型
prob = cp.Problem(objective, constraints)
prob.solve(solver=cp.CBC, verbose=False, maximumSeconds=100)
print("status:", prob.status)
print("最优x值:", x.value)
print("S_b:", S_b.value, "S_c:", S_c.value, "S_d:", S_d.value)

约束3的可行性分析 #

该约束完全可实现：

• 通过引入二进制辅助变量，我们将“判断最小值”的非线性逻辑转化为了线性约束，符合MILP的求解要求；
• 大M值取sum的最大可能值，确保约束在非激活状态下不会限制可行域；
• CVXPY的CBC求解器支持处理包含二进制变量的线性约束，能够正确求解该模型。

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