Question: 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!

Answer: here are some snippets that I would recommend to use in the model:  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()   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)   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)   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)   Don't forget to change model type to nonlinear:  model.params.FuncNonlinear = 1 model.update()   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.