在 Gurobi 中处理对数和指数目标函数

Question

我正在处理一个非常复杂的优化问题，该问题基本上使用医学参数来计算某个人患心血管疾病的风险。当然，我想最小化这种概率，这是由以下目标函数描述的：

地点：

和 beta 将始终是浮点数。

x向量也是一个浮点数，并且有一个初始值，代表该人当前的医疗状况。

现在我的问题是，考虑到对数和指数的存在，如何将目标函数正确加载到我的 Gurobi 模型中？以下是我当前的代码：

def calculate_beta(x2, x3, i):
    beta = np.zeros((2, 2, 15))
    beta[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]
    beta[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]
    beta[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]
    beta[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]
    return beta[x2, x3, (i-1)]

def calculate_xi(x1, x2, x3, x5, x6, x8, x9, x10):
    xi = ((np.log(x1) * calculate_beta(x2, x3, 1))
          + ((np.log(x1)**2) * calculate_beta(x2, x3, 2))
          + (np.log(x9) * calculate_beta(x2, x3, 3))
          + (np.log(x1) * np.log(x9) * calculate_beta(x2, x3, 4))
          + (np.log(x8) * calculate_beta(x2, x3, 5))
          + (np.log(x1) * np.log(x8) * calculate_beta(x2, x3, 6))
          + (np.log(x10) * calculate_beta(x2, x3, 7))
          + (np.log(x1) * np.log(x10) * calculate_beta(x2, x3, 8))
          + (0.0 * calculate_beta(x2, x3, 9))
          + (0.0 * calculate_beta(x2, x3, 10))
          + (x6 * calculate_beta(x2, x3, 11))
          + (np.log(x1) * x6 * calculate_beta(x2, x3, 12))
          + (x5 * calculate_beta(x2, x3, 13)))
    return xi

def ascvd_risk(x1, x2, x3, x5, x6, x8, x9, x10):
    probability = 1 - calculate_beta(x2, x3, 15) ** (np.exp(calculate_xi(x1, x2, x3, x5, x6, x8, x9, x10) - calculate_beta(x2, x3, 14)))
    return probability

model = Model("ascvd_risk")

x1 = model.addVar(vtype=GRB.INTEGER, name="x_1_age")
x2 = model.addVar(vtype=GRB.BINARY, name="x_2_gender")
...
x23 = model.addVar(vtype=GRB.BINARY, name="x_23_daily_alcohol")
x24 = model.addVar(vtype=GRB.CONTINUOUS, name="x_24_bmd")

model.setObjective(ascvd_risk(x1, x2, x3, x5, x6, x8, x9, x10), GRB.MINIMIZE)

现在，当我调用最后一行代码时，出现以下错误：

IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices

可以追溯到

----> 9 return beta[x2, x3, (i-1)]

从我在 Gurobi 论坛中看到的内容来看，我必须将其近似为分段线性函数，但我不知道具体该怎么做。

提前致谢！

Answer 1

以下是我建议在模型中使用的一些片段：

1. 将年龄变量转换为二进制：

from gurobipy import Model, GRB, Var, quicksum

# age variables, ages from 20 to 40
x1 = {}
for i in range(20,40):
    var = model.addVar(vtype=GRB.BINARY, name=f"x_1_age_{i}")
    model.update()
    x1[var] = i
    
# With SOS constraint you will only have single non-zero age_var
model.addSOS(GRB.SOS_TYPE1, [x_age for x_age in x1.keys()])
model.update()

2. 计算对数而不添加新的非线性

from gurobipy import Model, GRB, Var, quicksum

def var_log(x_age: dict[Var, int]):
    result = [x_var * np.log(age) for x_var, age in x_age.items()]
    return quicksum(result)

tt = var_log(x1)

3. 从数组中计算 beta：

beta = np.zeros((2, 2, 15))
beta[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]
beta[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]
beta[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]
beta[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]
beta = np.repeat(beta, 20, axis=1)


def calculate_beta(x_gender: Var, x_age: dict[Var, int], j):
    global beta

    beta_male = quicksum(x_var * beta[0,age-20, j] for x_var, age in x_age.items())
    beta_female = quicksum(x_var * beta[1,age-20, j] for x_var, age in x_age.items())

    return x_gender * beta_male + (1- x_gender) * beta_female

tt = calculate_beta(x2, x1, 5)