Handle logarithmic and exponential objective function in Gurobi

I am working with a very complex optimization problem which basically uses medical parameters to calculate the risk of a certain person developing a cardiovascular disease. I want to, of course, minimize such probability, which is described by the following objective function:

where:

and beta will always be a float.

The x vector is also a float and will have an initial value, which represents the person's current medical conditions.

Now my question is, how can I load my objective function correctly into my Gurobi Model, given the presence of logarithms and exponentials? Below is my current code:

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)

Now, when I call the last line of code, I get the following error:

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

which traces back to

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

From what I've seen in the Gurobi forums, I have to approximate it to a piece-wise linear function, but I don't know how exactly I can do it.

Thanks in advance!

here are some snippets that I would recommend to use in the model:

1. Split variable X3 into multiple binary variables. This is need to later get variable specific constants from the beta matrix:

from gurobipy import Model, GRB, Var, quicksum

# race variables...
# Map: gurobi_var -> var_index(int)
x3 = {}
for i in range(20):
    var = model.addVar(vtype=GRB.BINARY, name=f"x_3_race_{i}")
    model.update()
    x3[var] = i
    
# With SOS constraint you will only have single non-zero X3
model.addSOS(GRB.SOS_TYPE1, [x for x in x3.keys()])
model.update()

2. Non-linearity must be implemented using Gurobi non-linear constraints:

from gurobipy import Model, GRB, Var, quicksum


def var_log(x_var: Var):
    y = model.addVar(vtype=GRB.CONTINUOUS, name=f"y_nonlinear")
    gc = model.addGenConstrLog(x_var, y)
    model.update()
    
    # inform Gurobi that constraint shouldn't be treated as piecewise-linear approximations
    gc.FuncNonlinear = 1
    return y

tt = var_log(x2)

3. This is an example how to boradcast beta constants into gurobi expressions:

# --> generate random data 
# you can skip this part

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)

# <-- end of random data generation


def calculate_beta(x_gender: Var, x_3: dict[Var, int], j) -> Var:
    global beta

    beta_male = quicksum(x_var * beta[0,index, j] for x_var, index in x_3.items())
    beta_female = quicksum(x_var * beta[1,index, j] for x_var, index in x_3.items())

    beta_full = x_gender * beta_male + (1- x_gender) * beta_female

    # return gurobi variable class instead of Gurobi expression

    y = model.addVar(vtype=GRB.CONTINUOUS)
    gc = model.addConstr(y == beta_full)
    model.update()

    return y

tt = calculate_beta(x2, x3, 5)

4. Implement ASVCD using Gurobi non-linear constraints.

Link redirects to documentation on such constraints: https://www.gurobi.com/documentation/current/refman/general_constraints.html#subsubsection:GenConstrFunction

# ---> Dummy function instead of Xi
# it is needed to check that ASVCD constraints will compile 
# Replace with your implementation

calculate_xi = lambda x1, x2, x3, x5, x6, x8, x9, x10: x2

# <--- end of dummy function


def ascvd_risk(x1, x2, x3, x5, x6, x8, x9, x10):
    """
        prob = 1 - exp(np.ln(beta) * (exp(xi) - beta))
    """

    def get_var_exp(x_var: Var) -> Var:

        """
            Creates new var that is equal to `exp(x_var)` (same as `2.7 ** x_var`)
        """

        y = model.addVar(vtype=GRB.CONTINUOUS, name=f"xi_var")
        gc = model.addGenConstrExp(x_var, y)    
        model.update()
        gc.FuncNonlinear = 1
        return y

    def get_var_log(x_var: Var) -> Var:
        """
            Creates new var that is equal to `log(x_var)`
        """

        y = model.addVar(vtype=GRB.CONTINUOUS)
        gc = model.addGenConstrLog(x_var, y)
        model.update()
        gc.FuncNonlinear = 1
        return y

    var_beta_log = get_var_log(calculate_beta(x2, x3, 14))
    var_xi = calculate_xi(x1, x2, x3, x5, x6, x8, x9, x10)
    
    var_power = model.addVar(vtype=GRB.CONTINUOUS)
    model.addConstr(var_power == var_beta_log * (calculate_beta(var_xi) - get_var_beta(x2, x3, 13)))
    probability = 1 - get_var_exp(var_power)
    return probability

prob = ascvd_risk(x1, x2, x3, x5, x6, x8, x9, x10)

5. Don't forget to change model type to nonlinear:

model.params.FuncNonlinear = 1
model.update()

6. Off topic and other thoughts. As an alternative to Gurobi (if you won't be able to pass all your constraints to the model) I would suggest to have a closer look on differential evolution algos like scipy.optimize.differential_evolution. They work very well with non-linear systems that have less than 40-50 variables.