2万数据点值优化问题：寻求SciPy线性规划自动化方案

问题背景 #

现有20000个数据点(Xi)，每个数据点包含一个正值及4个属性：

• 属性A：99个互斥类别，每个数据点仅属于其中一类
• 属性B、C、D：支持将数据点的值按固定比例分配至多个类别
  已知各数据点的历史值，需求解New Value列，满足以下约束与目标。

硬约束 #

• New Value必须为正（数值计算中可替换为≥极小值ε，如1e-6）
• 所有数据点的New Value之和等于已知总值
• 各属性下，所有数据点对应类别分配比例与New Value的乘积之和，等于该类别目标值；且每个属性下所有类别目标值之和等于已知总值（共152个目标值，部分为关键约束）

软约束与优化目标 #

• 优先让Diff = New Value - Prev Value为正
• 小值数据点允许更大幅度的变动
• 核心目标：最小化每行Diff的相对变化幅度，均匀分配调整比例，避免单点大幅变动

现有尝试 #

当前每年通过数千次手动迭代解决，已尝试Python自动化及Matlab fmincon，现寻求同类问题资源，以及基于SciPy.optimize.linprog/milp的线性规划建模方案与优化建议。

基于SciPy的建模方案 #

一、约束线性化整理 #

变量定义

设x_i为第i个数据点的New Value，作为决策变量（共20000个）。

硬约束转化

1. 非负约束：x_i ≥ 1e-6（替代严格正，满足数值计算要求）
2. 总值约束：Σ(x_i) = Total（已知总值，若属性目标值之和已等于Total可省略）
3. 属性A约束：对每个A类k，Σ(x_i for i ∈ A_k) = Target_Ak（A_k为该类别包含的数据点索引集合）
4. 属性B/C/D约束：对每个B类m，Σ(x_i * p_ibm for all i) = Target_Bm（p_ibm为数据点i分配至B类m的固定比例，同理处理C/D类）

软约束转化

• 优先Diff为正：可设x_i ≥ Prev_i作为硬约束（若无解则放宽），或引入惩罚项优先满足该条件
• 小值数据点更自由：给每个数据点的变化量添加权重，权重与历史值正相关（如w_i = Prev_i），降低小值点变化对目标的影响

二、线性规划（linprog）实现：最小化加权绝对偏差 #

通过引入辅助变量将绝对值目标转化为线性约束，适合用linprog求解：

import numpy as np
from scipy.optimize import linprog

# 已知数据初始化
n = 20000
prev_values = np.array([...])  # 历史值数组
total = ...  # 已知总值
# 属性A约束构造（示例）
A_eq = []
b_eq = []
for k in range(99):
    row = np.zeros(n)
    row[category_A_indices[k]] = 1  # category_A_indices[k]为A类k的索引列表
    A_eq.append(row)
    b_eq.append(target_A[k])
# 属性B/C/D约束构造（示例）
for m in range(len(target_B)):
    row = np.zeros(n)
    for i in range(n):
        row[i] = p_ibm[i][m]  # p_ibm为数据点i到B类m的比例矩阵
    A_eq.append(row)
    b_eq.append(target_B[m])
# 合并约束为数组
A_eq = np.array(A_eq)
b_eq = np.array(b_eq)

# 目标：加权绝对偏差，小值点权重低
# 变量扩展为[x_1,...,x_n, d_1,...,d_n]，d_i为绝对偏差辅助变量
c = np.hstack([np.zeros(n), prev_values])  # d_i的权重与历史值正相关

# 绝对偏差约束：d_i ≥ |x_i - prev_values[i]|
A_ub = np.vstack([
    np.hstack([-np.eye(n), np.eye(n)]),  # -x_i + d_i ≥ -prev_values[i]
    np.hstack([np.eye(n), np.eye(n)])    # x_i + d_i ≥ prev_values[i]
])
b_ub = np.hstack([-prev_values, prev_values])

# 非负约束：x_i≥1e-6，d_i≥0
bounds = [(1e-6, None)] * n + [(0, None)] * n

# 扩展等式约束到辅助变量
A_eq_lp = np.hstack([A_eq, np.zeros((A_eq.shape[0], n))])
b_eq_lp = b_eq

# 求解
result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq_lp, b_eq=b_eq_lp, bounds=bounds, method='highs')
new_values = result.x[:n]

三、混合整数线性规划（milp）实现：优先满足Diff为正 #

引入0-1变量标记未满足Diff>0的点，通过惩罚项优先满足该软约束：

from scipy.optimize import milp, LinearConstraint, Bounds

# 变量：x_i(连续), d_i(连续), y_i(0/1整数)
n_vars = 3 * n
# 目标：加权绝对偏差 + 未满足Diff>0的惩罚（惩罚系数P=1e5）
c = np.hstack([np.zeros(n), prev_values, np.full(n, 1e5)])

# 约束构造
constraints = []
# 绝对偏差约束
constraints.append(LinearConstraint(
    np.hstack([-np.eye(n), np.eye(n), np.zeros((n, n))]),
    -prev_values, np.inf
))
constraints.append(LinearConstraint(
    np.hstack([np.eye(n), np.eye(n), np.zeros((n, n))]),
    prev_values, np.inf
))
# Diff>0软约束：x_i ≥ prev_i - M*y_i（M为足够大的数，如total）
M = total
constraints.append(LinearConstraint(
    np.hstack([np.eye(n), np.zeros((n, n)), -M*np.eye(n)]),
    prev_values, np.inf
))
# 硬约束：属性目标匹配
constraints.append(LinearConstraint(
    np.hstack([A_eq, np.zeros((A_eq.shape[0], 2*n))]),
    b_eq, b_eq
))

# 变量边界
bounds = Bounds(
    lb=np.hstack([np.full(n, 1e-6), np.zeros(n), np.zeros(n)]),
    ub=np.hstack([np.full(n, np.inf), np.full(n, np.inf), np.full(n, 1)])
)
# 标记整数变量（y_i为0/1）
integrality = np.hstack([np.zeros(n), np.zeros(n), np.ones(n)])

# 求解
result = milp(c, constraints=constraints, bounds=bounds, integrality=integrality)
new_values = result.x[:n]

同类问题资源与优化建议 #

1. 问题归类：属于多维度约束下的资源分配/数值校准问题，常见于财务预算拆解、销售目标分配、供应链需求校准等场景。
2. 优化建议：
   • 约束预处理：先验证各属性目标值之和是否等于总值，若不匹配需先调整目标值（如按比例缩放非关键约束），否则问题无解。
   • 初始值优化：用历史值作为初始解，计算约束偏差后优先调整小值数据点，减少优化迭代成本。
   • 非线性目标替代：若需更精准的均匀百分比变化控制，可改用二次规划（如scipy.optimize.minimize的SLSQP方法，目标设为Σ((x_i/prev_i - 1)^2)），或用CVXPY建模二次规划，灵活性更强。
   • 规模拆分：将问题按属性A的类别拆分，先分配总值到各A类，再在类内处理B/C/D约束，降低单模型变量规模。
3. 参考方向：可参考运筹学教材中“运输问题”“分配问题”的扩展应用，以及SciPy官方文档中linprog/milp的多约束案例。

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