帮助格式化和解决一个非常大的方程组问题
这是一个现实中的问题，以前通过成千上万次的手动迭代解决过。
请原谅我对StackOverflow的不熟悉以及提问格式的不足。Dave在评论中正确地重写了我的问题，我将在下面引用他的话。如果我能找到方法，我也想感谢他。
如果你能让问题更集中，先定义术语再使用，并且明确输入是什么，可能会得到更好的回应。例如，“一个数据点是一个容器，里面有一个叫做‘值’的正整数，一个叫做‘A’的属性，里面有一个整数，还有属性B、C、D，每个属性里面都有整数和浮点数的组合，浮点数的总和为1.0。我的输入大约有2万个这样的数据点，我的目标是找到这些点的新值，保持其他不变，最大化（新值 - 旧值）。”
我有大约20,000个数据点（Xi）。每个数据点都有一个大于零的值。每个数据点还有4个属性。属性A有99个可能的类别，每个数据点只能属于其中一个。剩下的三个属性可以将值分配到多个类别。例如：Xi的80%值属于属性B的类别2，而剩下的20%属于属性B的类别5。
我还有每个数据点在过去几年的值（每年一个值）。
| 点 | 之前的值 | 新值 | 差异 | 属性A | 属性B | 属性C | 属性D |
|---|---|---|---|---|---|---|---|
| X1 | 68 | 72 | 4 | 1: 100% | 2: 80% | 3: 100% | 7: 90% |
| 5: 20% | 9: 10% | ||||||
| X2 | 56,000 | 66,000 | 10,000 | 7: 100% | 1: 50% | 3: 90% | 2: 100% |
| 5: 50% | 6: 10% |
我需要解决新值这一列，需遵循以下约束条件：
新值必须是正数
所有(Xi)的总和 = 已知值
所有(Xi) * 每个属性分配给类别的百分比 = 目标值
所有目标值的总和 = 已知值
属性A
A1: 目标值
A2: 目标值
...
A99: 目标值
A1 + A2 +... A99 = 已知值
属性B
B1: 目标值
B2: 目标值
...
B27: 目标值
B1 + B2 +... B27 = 已知值
属性C
C1: 目标值
C2: 目标值
...
C18: 目标值
C1 + C2 +... C18 = 已知值
属性D
D1: 目标值
D2: 目标值
...
D8: 目标值
D1 + D2 +... D8 = 已知值
软约束
差异列中的值应该是正数
较小的值应该有更多的变化自由度，而较大的值则应限制变化
最终目标
尽量减少每行的差异列的值。在查看与之前值的百分比变化时，尽量均匀分配变化，而不是全部归因于单一行。
如前所述，这个问题每年通过手动迭代解决一次。我们开发了一种工作流程，首先通过将之前的值转换为总值的百分比来填充起始值，然后将其乘以新的已知值，这个值总是大于之前的总值。我们使用BI工具同时查看起始值和目标值之间的所有总差异。然后，手动迭代以[从X34减去10,000并加到X452]的形式添加到一个运行列表中，并更新BI工具以显示新状态。这个运行列表可能有成千上万行。此外，我所说的“解决”是指我们最终可以生成一个符合所有约束的解决方案，但我们清楚这只是一个解决方案，而不是最佳解决方案。
我们已经尝试通过Python自动化迭代过程，取得了一定的成功。我们还与Matlab的代表进行了交谈，他们相信可以使用fmincon接近解决方案，我们仍然可能会追求这个方向，但我想探索其他选择。
我在这篇文章中请求的并不一定是一个解决方案（虽然我会接受），而是类似问题的资源链接（一个元素的变化导致其他元素的变化）。或者可能是帮助我更清晰地定义数学问题。我已经研究了几种优化和遗传算法，但似乎没有一个符合我的需求。
- “所有(Xi)的总和 * 每个属性分配给类别的百分比 = 目标值”这个问题是一个整体的目标值，还是每个类别的特定目标值，或者其他什么？
每个属性的每个类别都有一个目标值。属性D有8个目标值。这8个值的总和等于属性C的18个目标值的总和（或任何属性）。这是一个封闭系统。总共有99 + 27 + 18 + 8 = 152个目标值。实际上，这些目标值中有些比其他的更重要。属性B必须匹配所有27个目标，而属性C中只有5个目标是关键的。
- 你优化的参数是什么？例如，你能改变类别之间的百分比分配吗？ – Nick ODell
实际上，百分比可以稍微变化。比如从80/20变成75/25。老实说，我们一直保持这些分配不变，以简化已经复杂的问题。你不能做的是添加一个不存在的类别。例如，将80/20改为70/20/10。我们希望优化（最小化）的参数是之前值和新值之间的差异。同时，理想的解决方案也应该是平滑的。给出两个解决方案：第一个解决方案中，大多数数据点变化很小，只有一个数据点变化很大。而在第二个解决方案中，所有数据点变化较大，但没有单个点变化很大。第二个解决方案更受欢迎。
如果我遗漏了上面提到的内容，请原谅，但到目前为止，所有约束似乎都是线性的。如果目标是最小化差异的绝对值之和，那么这可以表示为线性规划问题，并使用scipy.optimize.linprog或scipy.optimize.milp解决，这比minimize更适合处理大问题。如果你对这种方法感兴趣，我可以提供一个示例，帮助你以类似的方式构建问题。 – Matt Haberland
@Matt 是的，我非常想看看你如何在SciPy中构建这个问题的解决方案。可能还会结合Cary Swoveland和Nick ODell的建议。
1 个回答
注意：我没有更多时间来润色这个内容。即使第一次读起来有点难懂，我希望它至少能让你了解是否想自己去自动化这个解决方案，或者想雇人来做。
线性规划的一个好处是，如果它足够灵活，可以解决你的问题，那么它的解决方案非常高效。基本上，你可以保证（除非遇到一些特殊的数值问题）得到一个全局最优解，或者得到一个不可行或无界的证明。
scipy.optimize.milp 可以解决线性和混合整数线性规划问题，它的界面更清晰，与底层求解器（HiGHS）的连接也更紧密，所以我建议使用这个。下面的所有代码都是用普通的 NumPy 数组写的，但你最好使用 稀疏数组，因为数组会很大，很多元素会是零。
我有大约 20,000 个数据点（Xi）。
在优化的术语中，这些数据点的新值将被称为“决策变量”。假设你有 mx 个这样的变量。
听起来你的“目标值”Aj、Bj、Cj 和 Dj 也是需要由求解器选择的变量，因为你用它们来写约束，所以这些也算是“决策变量”。假设你分别有 ma=99、mb=27、mc=18 和 md=8。
所以现在，假设你有 m = mx + ma + mb + mc + md 个决策变量，我们把这些决策变量统称为 x。
新值必须是正数
默认情况下，所有决策变量都被限制为 非负。没有办法指定严格的不等式约束（“正数”），但你可以用一个 Bounds 对象来表示包含边界，以确保决策变量满足 lb <= x <= ub。
# vector of size `m` of lower bounds, one for each decision variable
# may be a small nonzero number to ensure that values are strictly positive
lb = ...
# vector of size `m` of upper bounds, one for each decision variable
# may be `np.inf` to indicate "no upper bound"
ub = ...
bounds = Bounds(lb, ub)
（另一种策略是希望所有值在最优解中都是非零的，如果发现某些值恰好为零，再添加边界。）
所有 (Xi) 的总和 = 已知值
我们将用一个 LinearConstraint 来指定这个约束，它表示的形式是 lb <= A @ x <= ub。
A1 = np.zeros((1, m))
A1[:mx] = 1 # you only want to sum your Xi here, not your Aj, Bj, etc...
b1 = np.asarray([...]) # your known value as a vector with one element
constraint1 = LinearConstraint(A1, lb=b1, ub=b1)
所有目标值的总和按属性分配的百分比 = 目标值
实际上，百分比可能会稍微变化。可能不是 80/20，而是 75/25。老实说，我们一直保持这些比例不变，以简化一个已经复杂的问题。
我假设这些百分比是已知的。如果不是，这些百分比将成为决策变量，你的约束将涉及决策变量的乘积，这将使你的问题变为非线性（在决策变量中）。非线性问题不能用 linprog 或 milp 来解决。例外的是，我们可以用二进制决策变量将数据点分配到 A 的单一类别。我会在最后讨论这个，但现在假设你已经完成了所有的分配。
这是我认为你可以如何以 lb <= A @ x <= ub 的形式指定这个约束。
# Matrix with 152 rows and `mx` columns specifying percentages in each category
# Assume for now that you will manually assign datapoints to categories for
# attribute A; we'll consider that later
A2a = ...
# To ensure that the sums of these Xi * percentages equals the target value,
# we subtract the target value...
A2b = -np.eye(ma + mb + mc + md)
A2 = np.concatenate([A2a, A2b], axis=1)
# ... and the result must be zero for each row.
b2 = np.zeros(ma + mb + mc + md)
constraint2 = LinearConstraint(A2, lb=b2, ub=b2)
每个类别的所有目标值的总和 = 已知值
希望这开始变得有意义：
A3 = np.zeros((4, m)) # one constraint per category
A3[0, mx:mx + ma] = 1 # sum all the Aj in the zeroth row
A3[1, mx + ma: mx + ma + mb] = 1 # sum all the Bj in the next row
...
b3 = ... # target values for each sum
较小的值应该有更多的自由度去变化，而不是较大的值
你可以重新表述这个问题，不再使用 Xi 本身作为决策变量，而是将决策变量设为 Xi 与旧值的 比率。假设你这样做了。
Diff 列中的值应该是正数
为了简单起见，我们将这个设为一个硬约束：每个 Xi 比率决策变量的下限为 1。（软约束可以融入目标函数，但这取决于你的创造力。）
寻求按行最小化 Diff 列中的值。
再次强调，决策变量现在是 Xi 的 比率，你将所有这些比率的下限设为 1。现在的目标是 最小化最大比率。为此，我们添加一个新的决策变量，要求它大于所有 Xi 比率的最大值。我们通过限制这个新决策变量大于 每个 Xi 比率来确保这一点：
A = np.zeros((mx, m + 1))
A[:, -1] = 1 # corresponding with the new decision variable
i = np.arange(mx)
A[i, i] = -1 # subtract ratio `i` in row `i`
lb = np.zeros(mx) # must be greater than zero
ub = np.full(mx, np.inf)
你的目标函数将变为：
c = np.zeros(m + 1)
c[-1] = 1 # minimize the new decision variable
注意：因为你添加了这个新的决策变量，所有约束矩阵都需要多一列。所有约束矩阵需要有相同数量的列，也就是决策变量的总数。
注意 2：我将“Diff 列中的值应该是正数”设为一个硬约束。如果这不是一个约束，你表达目标函数的方式将需要改变；它不能再是“最小化最大 Xi 比率”，因为某些比率将小于 1。可以推测你也希望确保它们不会 远远 小于 1。你可能希望目标函数变为“最小化 比率与 1 之间差值的绝对值的最大值”。我不会在这里尝试描述如何适应这一点，但可以参考 AIMMS建模指南中的“绝对值”和“最小最大目标”技巧。
现在我们回到数据点分配到属性 A 的类别。基本思路是这样的。你添加 mx * ma 个 新决策变量（是的，像 2,000,000 个），每个变量要么是 1（如果 Xi 被分配到类别 Aj），要么是 0（否则）。你可以通过 milp 的 integrality 参数确保这些决策变量是“二进制”的。你添加的线性约束与上面类似，形成一个矩阵 A，它有 mx 行和 m + 1 + mx * ma 列。这个矩阵几乎每个元素都是 0，但在每一行 j 中，会有 ma=99 个 1，对应于将数据点 i 分配到类别 j 的二进制决策变量。你将所有 lb 和 ub 的元素都设为 1——与单个数据点对应的所有二进制决策变量的总和必须等于 1。（这确保了一个决策变量只能属于一个类别。）
在 SciPy 用户指南 中有一些线性规划（LP）和混合整数线性规划（MILP）的教程，你还可以在 Medium 和其他地方找到其他教程。