如何在Scipy.optimize.minimize中设置非零变量数量限制约束？

你的核心需求是在投资组合优化中限制非零权重股票数量不超过400，这属于基数约束（Cardinality Constraint），这类约束是非线性且非凸的，Scipy.optimize.minimize的SLSQP这类连续优化方法原生无法处理，必须将问题建模为混合整数二次规划（MIQP）——因为你的目标函数是平方和（二次），同时引入二进制整数变量来标记非零权重。

1. 问题建模转换 #

首先把原问题转化为MIQP形式：

变量定义 #

• 连续变量 w_i：第i只股票的投资组合权重，满足 0 ≤ w_i ≤ 1
• 二进制变量 z_i ∈ {0,1}：标记第i只股票是否取非零权重，z_i=1 表示 w_i>0，z_i=0 表示 w_i=0

约束条件 #

• 基数约束：限制非零权重的股票数量 ≤400

  ∑z_i ≤ 400 （i从1到700）
  

• 关联约束：保证二进制变量和权重变量的逻辑一致性

  w_i ≤ z_i （因为w_i的上界是1，所以用1作为缩放系数，确保z_i=0时w_i必须为0；z_i=1时w_i可在0-1范围内取值）
  

  可选补充（避免数值误差导致的极小非零值）：

  w_i ≥ ε*z_i （ε取极小值，比如1e-6，确保z_i=1时w_i有实际意义的非零值）
  

• 你已实现的行业分配等其他约束（直接保留原有形式即可）

目标函数 #

保留你原有的跟踪误差最小化目标：

minimize ∑(w_i - b_i)² （b_i为第i只股票的基准指数权重）

2. 工具选择与代码示例 #

Scipy不支持整数变量，需要用专门的混合整数规划建模库，推荐使用PuLP（易上手，开源免费），或商业求解器Gurobi/Cplex（性能更优，适合大规模问题）。以下是PuLP的实现示例：

import pulp

# 模拟700只股票的场景
stock_count = 700
max_non_zero_stocks = 400
# 替换为你的实际基准权重数据
benchmark_weights = [0.001 + 0.0001*i for i in range(stock_count)]

# 创建MIQP问题
prob = pulp.LpProblem("Portfolio_Tracking_Optimization", pulp.LpMinimize)

# 定义变量：投资组合权重（连续变量）
weights = pulp.LpVariable.dicts(
    "StockWeight", 
    range(stock_count), 
    lowBound=0, 
    upBound=1, 
    cat="Continuous"
)
# 定义二进制变量：标记是否为非零权重
indicators = pulp.LpVariable.dicts(
    "NonZeroIndicator", 
    range(stock_count), 
    cat="Binary"
)

# 目标函数：最小化跟踪误差平方和
prob += pulp.lpSum(
    [(weights[i] - benchmark_weights[i])**2 for i in range(stock_count)]
)

# 核心基数约束：非零权重股票数量≤400
prob += pulp.lpSum([indicators[i] for i in range(stock_count)]) <= max_non_zero_stocks

# 关联权重与标记变量的约束
eps = 1e-6  # 避免数值误差的极小值
for i in range(stock_count):
    prob += weights[i] <= indicators[i]
    prob += weights[i] >= eps * indicators[i]

# 添加你已实现的行业分配约束（示例）
# 比如：行业A包含股票0-149，权重和需等于基准行业权重0.2
# industry_a_stocks = range(150)
# prob += pulp.lpSum([weights[i] for i in industry_a_stocks]) == 0.2

# 求解问题（使用CBC开源求解器，支持MIQP）
solver = pulp.PULP_CBC_CMD(msg=1)  # msg=1显示求解过程日志
prob.solve(solver)

# 输出结果统计
non_zero_count = 0
for i in range(stock_count):
    weight_val = pulp.value(weights[i])
    if weight_val > eps:
        non_zero_count += 1
        # 可选：打印每只非零权重股票的数值
        # print(f"股票{i}权重: {weight_val:.6f}")
print(f"最终非零权重股票数量: {non_zero_count}")
print(f"目标函数值（跟踪误差平方和）: {pulp.value(prob.objective):.6f}")

3. 关键注意事项 #

• 数值稳定性：添加w_i ≥ ε*z_i可以避免求解器因数值精度问题生成极小的非零权重，ε可根据实际需求调整
• 求解性能：700个连续变量+700个二进制变量的规模，开源CBC求解器足够处理；如果求解速度慢，可尝试商业求解器Gurobi/Cplex，它们对MIQP的优化更高效
• 约束优先级：如果行业分配等约束与基数约束冲突，求解器会返回不可行解，需要调整约束的合理性

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