如何在Scipy中实现食品营养组合优化的线性规划模型

我有多组食品类别A、B……，每组食品的营养信息以向量x_a1、x_a2……、x_b1、x_b2……表示；另有目标营养值向量y。希望从这些食品集中选择组合，使目标营养y与所选食品营养向量之和的差异最小化：
$$\min_x ||x−y||2$$
其中x为所有所选食品的元素-wise求和：
$$x = \sum_i(x{ij})$$
$i∈{A,B,C,...}$，$j∈{1,2,3,...}$。

为通过线性规划解决该问题，建模如下：设Z为m行n列的矩阵（n种食品，m种营养成分），w_j为二进制向量。通过引入松弛变量t_i和非负剩余变量s_i，将最小化L1范数：
$$\sum_i | \sum_j(Z_{ij}w_j) - y |$$
转化为线性问题：
$$\min \sum_i(s_i + t_i)$$
约束条件：
$$\sum_j(Z_{ij}w_j) - s_i + t_i = y \quad \text{for all } i$$
$$\sum_j(w_j) >= 1 \quad \text{for } j∈{A,B,C,...}$$

以下是在Scipy中实现该线性规划模型的步骤：

Scipy的scipy.optimize.linprog工具可解决线性规划问题，注意linprog默认支持最小化问题，且变量非负；由于w_j是二进制变量，需使用支持整数约束的method='highs'求解器。

1. 构造目标函数系数向量c #

目标是最小化$\sum_i(s_i + t_i)$，总变量包含：

• n个二进制变量w_j（对应n种食品）
• m个非负剩余变量s_i
• m个非负松弛变量t_i

目标函数中，w_j的系数为0，s_i和t_i的系数为1，示例构造代码：

import numpy as np
from scipy.optimize import linprog

# 自定义参数：m种营养成分，n种食品
m = 3  # 示例：蛋白质、碳水、脂肪
n = 5  # 示例：5种可选食品

# 构造目标系数向量c：前n位对应w_j，中间m位对应s_i，最后m位对应t_i
c = np.concatenate([np.zeros(n), np.ones(m), np.ones(m)])

2. 构造约束条件矩阵与向量 #

约束分为营养平衡约束和选品数量约束两类：

营养平衡约束 #

对每种营养i，约束式$\sum_j Z_{ij}w_j - s_i + t_i = y_i$需转化为等式约束矩阵A_eq和向量b_eq：

# 示例营养矩阵Z和目标向量y，实际替换为你的真实数据
Z = np.random.rand(m, n)  # m行n列，每行对应一种营养的各食品含量
y = np.random.rand(m)     # m维目标营养向量

# 构造等式约束矩阵A_eq
A_eq = np.zeros((m, n + 2*m))
for i in range(m):
    A_eq[i, :n] = Z[i, :]       # 对应w_j的系数
    A_eq[i, n + i] = -1         # 对应s_i的系数
    A_eq[i, n + m + i] = 1      # 对应t_i的系数

b_eq = y.copy()  # 等式约束的右侧向量

至少选一种食品的约束 #

约束式$\sum_j w_j >= 1$需转换为linprog支持的A_ub @ x <= b_ub形式：

# 构造不等式约束：sum(w_j) >=1 → -sum(w_j) <= -1
A_ub = np.zeros((1, n + 2*m))
A_ub[0, :n] = -1
b_ub = np.array([-1])

3. 设置变量边界与整数约束 #

• w_j是二进制变量，边界为(0, 1)且需标记为整数
• s_i和t_i是非负变量，边界为(0, None)

# 定义每个变量的上下界
bounds = []
# w_j的边界：二进制
for _ in range(n):
    bounds.append((0, 1))
# s_i的边界：非负
for _ in range(m):
    bounds.append((0, None))
# t_i的边界：非负
for _ in range(m):
    bounds.append((0, None))

# 标记整数变量：前n个变量（w_j）为整数，其余为连续变量
integrality = np.concatenate([np.ones(n), np.zeros(2*m)])

4. 调用求解器并解析结果 #

使用highs求解器执行整数线性规划：

# 求解线性规划
result = linprog(
    c=c,
    A_ub=A_ub,
    b_ub=b_ub,
    A_eq=A_eq,
    b_eq=b_eq,
    bounds=bounds,
    integrality=integrality,
    method='highs'
)

# 解析输出结果
print("求解状态：", result.message)
print("最小L1差异值：", result.fun)
# 提取选中的食品（二进制变量接近1即为选中）
selected_food_indices = np.where(result.x[:n] >= 0.99)[0]
print("选中的食品索引：", selected_food_indices)
# 计算选中食品的营养总和
selected_nutrition = Z.T @ result.x[:n]
print("实际营养总和：", selected_nutrition)
print("目标营养值：", y)

注意事项 #

• 若你最初目标是最小化L2范数，当前建模转换为了L1范数的线性规划，二者优化目标不同。如需严格实现L2范数最小化，需使用二次规划工具（如scipy.optimize.minimize）。
• 确保营养矩阵Z与目标向量y的单位统一，避免数值差异过大导致求解不稳定。
• 当食品数量较多时，整数规划求解时间会增加，可通过调整options参数（如maxiter）优化求解效率。

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