在Gurobi中处理含对数与指数的目标函数问题

问题概述 #

正在处理一个心血管疾病（CV）风险最小化的优化问题，目标函数包含对数、指数及幂运算，核心结构如下：

• 风险概率公式：$P = 1 - \beta_{15}^{(x_2,x_3)} \exp\left( \xi - \beta_{14}^{(x_2,x_3)} \right)$
• 其中$\xi$是多变量加权组合：
  $$\begin{align*}
  \xi =& \beta_1 \ln x_1 + \beta_2 (\ln x_1)^2 + \beta_3 \ln x_9 + \beta_4 \ln x_1 \ln x_9 \
  &+ \beta_5 \ln x_8 + \beta_6 \ln x_1 \ln x_8 + \beta_7 \ln x_{10} + \beta_8 \ln x_1 \ln x_{10} \
  &+ \beta_{11} x_6 + \beta_{12} x_6 \ln x_1 + \beta_{13} x_5
  \end{align*}$$
  $\beta$是按二元变量$x_2$（性别）、$x_3$分组的预定义系数矩阵。

建模时触发IndexError，错误出在calculate_beta函数的索引行——因为$x_2$、$x_3$是Gurobi变量，无法直接作为NumPy数组的整数索引；同时目标函数中的非线性项（对数、指数、乘积）需用分段线性（PWL）近似处理。

第一步：修复索引错误 #

错误原因 #

$x_2$、$x_3$是Gurobi的Var对象（二元变量），不能直接用于NumPy数组索引，需通过逻辑约束+辅助变量实现系数组的选择。

实现方案 #

1. 预定义所有$\beta$系数，按$(x_2,x_3)$的四种组合存储：

import gurobipy as gp
from gurobipy import GRB
import math

beta_coeffs = {
    (0,0): [-29.799, 4.884, 13.540, -3.114, -13.578, 3.149, 2.019, 0.0, 1.957, 0.0, 7.574, -1.665, 0.661, -29.18, 0.9665],
    (0,1): [17.114, 0.0, 0.940, 0.0, -18.920, 4.475, 29.291, -6.432, 27.820, -6.087, 0.691, 0.0, 0.874, 86.61, 0.9533],
    (1,0): [12.344, 0.0, 11.853, -2.664, -7.990, 1.769, 1.797, 0.0, 1.764, 0.0, 7.837, -1.795, 0.658, 61.18, 0.9144],
    (1,1): [2.469, 0.0, 0.302, 0.0, -0.307, 0.0, 1.916, 0.0, 1.809, 0.0, 0.549, 0.0, 0.645, 19.54, 0.8954]
}

2. 引入二元辅助变量指示系数组选择，并添加逻辑约束：

model = gp.Model("ascvd_risk")
model.setParam('NonConvex', 2)  # 开启非凸二次项支持

# 原始变量定义
x1 = model.addVar(vtype=GRB.INTEGER, name="x_1_age", lb=20, ub=80)
x2 = model.addVar(vtype=GRB.BINARY, name="x_2_gender")
x3 = model.addVar(vtype=GRB.BINARY, name="x_3_group")
x5 = model.addVar(vtype=GRB.CONTINUOUS, name="x_5_hypertension", lb=0, ub=1)
x6 = model.addVar(vtype=GRB.CONTINUOUS, name="x_6_diabetes", lb=0, ub=1)
x8 = model.addVar(vtype=GRB.CONTINUOUS, name="x_8_cholesterol", lb=100, ub=300)
x9 = model.addVar(vtype=GRB.CONTINUOUS, name="x_9_bp", lb=80, ub=200)
x10 = model.addVar(vtype=GRB.CONTINUOUS, name="x_10_smoking", lb=0, ub=1)

# 系数组选择变量
z00 = model.addVar(vtype=GRB.BINARY, name="z_00")
z01 = model.addVar(vtype=GRB.BINARY, name="z_01")
z10 = model.addVar(vtype=GRB.BINARY, name="z_10")
z11 = model.addVar(vtype=GRB.BINARY, name="z_11")

# 互斥约束：仅能选择一组系数
model.addConstr(z00 + z01 + z10 + z11 == 1)

# 绑定z变量与x2、x3的逻辑关系
model.addConstr(z00 <= 1 - x2)
model.addConstr(z00 <= 1 - x3)
model.addConstr(z00 >= 1 - x2 - x3)

model.addConstr(z01 <= 1 - x2)
model.addConstr(z01 <= x3)
model.addConstr(z01 >= x3 - x2)

model.addConstr(z10 <= x2)
model.addConstr(z10 <= 1 - x3)
model.addConstr(z10 >= x2 - x3)

model.addConstr(z11 <= x2)
model.addConstr(z11 <= x3)
model.addConstr(z11 >= x2 + x3 - 1)

# 创建beta系数辅助变量
beta_vars = []
for k in range(15):
    b_var = model.addVar(vtype=GRB.CONTINUOUS, name=f"beta_{k+1}")
    beta_vars.append(b_var)
    model.addConstr(
        b_var == z00*beta_coeffs[(0,0)][k] + z01*beta_coeffs[(0,1)][k] + z10*beta_coeffs[(1,0)][k] + z11*beta_coeffs[(1,1)][k]
    )

第二步：分段线性近似处理非线性项 #

1. 对数函数$\ln(x)$ #

以$x_1$为例，生成采样点后用addGenConstrPWL创建近似约束：

# 近似ln(x1)
ln_x1 = model.addVar(vtype=GRB.CONTINUOUS, name="ln_x1")
x1_points = list(range(20, 81, 2))
ln_x1_points = [math.log(x) for x in x1_points]
model.addGenConstrPWL(x1, ln_x1, x1_points, ln_x1_points, name="pwl_ln_x1")

# 同理处理其他对数项
ln_x8 = model.addVar(vtype=GRB.CONTINUOUS, name="ln_x8")
x8_points = list(range(100, 301, 10))
ln_x8_points = [math.log(x) for x in x8_points]
model.addGenConstrPWL(x8, ln_x8, x8_points, ln_x8_points, name="pwl_ln_x8")

ln_x9 = model.addVar(vtype=GRB.CONTINUOUS, name="ln_x9")
x9_points = list(range(80, 201, 5))
ln_x9_points = [math.log(x) for x in x9_points]
model.addGenConstrPWL(x9, ln_x9, x9_points, ln_x9_points, name="pwl_ln_x9")

ln_x10 = model.addVar(vtype=GRB.CONTINUOUS, name="ln_x10")
# x10是0-1变量，单独处理采样点
x10_points = [0.01, 0.2, 0.4, 0.6, 0.8, 1]
ln_x10_points = [math.log(x) for x in x10_points]
model.addGenConstrPWL(x10, ln_x10, x10_points, ln_x10_points, name="pwl_ln_x10")

2. 平方项$(\ln x_1)^2$ #

基于$\ln x_1$的近似结果再次做分段线性近似：

ln_x1_sq = model.addVar(vtype=GRB.CONTINUOUS, name="ln_x1_sq")
y_points = [math.log(x) for x in range(20, 81, 2)]
y_sq_points = [y**2 for y in y_points]
model.addGenConstrPWL(ln_x1, ln_x1_sq, y_points, y_sq_points, name="pwl_ln_x1_sq")

3. 乘积项（如$\ln x_1 \times \ln x_9$） #

利用Gurobi的二次约束直接处理：

lnx1_lnx9 = model.addVar(vtype=GRB.CONTINUOUS, name="lnx1_lnx9")
model.addConstr(lnx1_lnx9 == ln_x1 * ln_x9, name="prod_lnx1_lnx9")

# 同理处理其他乘积项
lnx1_lnx8 = model.addVar(vtype=GRB.CONTINUOUS, name="lnx1_lnx8")
model.addConstr(lnx1_lnx8 == ln_x1 * ln_x8, name="prod_lnx1_lnx8")

lnx1_lnx10 = model.addVar(vtype=GRB.CONTINUOUS, name="lnx1_lnx10")
model.addConstr(lnx1_lnx10 == ln_x1 * ln_x10, name="prod_lnx1_lnx10")

lnx1_x6 = model.addVar(vtype=GRB.CONTINUOUS, name="lnx1_x6")
model.addConstr(lnx1_x6 == ln_x1 * x6, name="prod_lnx1_x6")

4. 指数函数$\exp(\xi - \beta_{14})$ #

先构建$\xi$的线性表达式，再对指数项做分段线性近似：

# 构建xi表达式
xi = (beta_vars[0] * ln_x1 + beta_vars[1] * ln_x1_sq + beta_vars[2] * ln_x9 + beta_vars[3] * lnx1_lnx9 +
      beta_vars[4] * ln_x8 + beta_vars[5] * lnx1_lnx8 + beta_vars[6] * ln_x10 + beta_vars[7] * lnx1_lnx10 +
      beta_vars[10] * x6 + beta_vars[11] * lnx1_x6 + beta_vars[12] * x5)

# 计算xi - beta_14
xi_minus_beta14 = model.addVar(vtype=GRB.CONTINUOUS, name="xi_minus_beta14")
model.addConstr(xi_minus_beta14 == xi - beta_vars[13], name="xi_minus_beta14")

# 近似exp(xi_minus_beta14)
exp_term = model.addVar(vtype=GRB.CONTINUOUS, name="exp_term")
z_points = [i for i in range(-50, 11)]
exp_z_points = [math.exp(z) for z in z_points]
model.addGenConstrPWL(xi_minus_beta14, exp_term, z_points, exp_z_points, name="pwl_exp")

5. 幂运算$\beta_{15}^{\exp(\cdot)}$ #

转化为$\exp(y \times \ln \beta_{15})$后处理：

# 近似ln(beta15)（beta15仅4个可能值，直接通过z变量绑定）
ln_beta15 = model.addVar(vtype=GRB.CONTINUOUS, name="ln_beta15")
model.addConstr(
    ln_beta15 == z00*math.log(0.9665) + z01*math.log(0.9533) + z10*math.log(0.9144) + z11*math.log(0.8954)
)

# 计算y * ln_beta15
y_ln_beta15 = model.addVar(vtype=GRB.CONTINUOUS, name="y_ln_beta15")
model.addConstr(y_ln_beta15 == exp_term * ln_beta15, name="prod_y_ln_beta15")

# 近似最终指数项
exp_final = model.addVar(vtype=GRB.CONTINUOUS, name="exp_final")
z2_points = [i for i in range(-2500, 1)]
exp_z2_points = [math.exp(z) for z in z2_points]
model.addGenConstrPWL(y_ln_beta15, exp_final, z2_points, exp_z2_points, name="pwl_exp_final")

第三步：设置目标函数并求解 #

# 目标是最小化1 - exp_final
model.setObjective(1 - exp_final, GRB.MINIMIZE)

model.optimize()

# 输出结果
if model.status == GRB.OPTIMAL:
    print("最优解：")
    for var in model.getVars():
        if var.x != 0:
            print(f"{var.varName}: {var.x:.4f}")
    print(f"最小风险概率：{model.objVal:.4f}")

关键注意事项 #

• 所有变量必须明确上下限，确保分段线性采样点覆盖实际取值范围。
• 采样点数量越多，近似精度越高，但模型求解速度会变慢，需根据需求权衡。
• 若不使用二次约束，乘积项也可通过分段线性近似实现，但会增加模型复杂度。

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