Search code examples
pythonvargurobinumerical-integration

How to iterate over a Gurobi decision variable in an integrated normal distribution that includes exponentiation


My problem: Iterating over Gurobi "Var" decision variable creates TypeError: '<' not supported between instances of 'Var' and 'int' and issue with exponentation (i.e. **/ pow())

Key information on the Gurobi optimization:

  • Objective function: for three items maximize the sum of (price * expected value)
  • Expected value is calculated through two defined formulas:
  • 1) PDF = probability density function
  • 2) EV = Expected value which is the integration of the PDF over a specific integral
  • The decision variable "upperBound" is supposed to maximize the upper bound of this integral, the lower bound is 0

Below the model:

from gurobipy import *
import pandas as pd, numpy as np, scipy.integrate as integrate
import math

mu = pd.DataFrame(np.array([10, 15]), index = ["product1", "product2"])
sigma = pd.DataFrame(np.array([1, 1.5]), index = mu.index.values)
price = pd.Series(np.array([10, 10]), index = mu.index.values)

m = Model("maxUpperBound")
ub = m.addVars(mu.index.values, vtype=GRB.INTEGER, name="upperBound")

def PDF (y, mu, sigma):
    return y * (1.0 / (sigma * (2.0 * math.pi)**(1/2)) * math.exp(-1.0 * (y - mu)**2 / (2.0 * (sigma**2))))

def EV(ub, mu, sigma):
    lb = 0
    value, error = integrate.quad(PDF, lb, ub, args=(mu, sigma))
    return value

m.setObjective(
    quicksum(
        price.loc[i] * EV(ub[i], mu.loc[i], sigma.loc[i])
        for i in mu.index.values
    ),
    GRB.MAXIMIZE)

m.addConstr(
    (quicksum(ub[i]
              for i in mu.index.values)
     <= 100),
    "Limit100"
)

m.optimize()

for v in m.getVars():
    print(v.varName, v.x)

Solution

  • Gurobi is not able to optimize with respect to the PDF function as this is not a mixed-integer problems comprised of linear, (convex or nonconvex) quadratic, or second-order constraints/objective.

    Instead solved this by calculating the expected value upfront for each product_ub combination. ub can takes values between 0 and 100 (see constraint). Used the expected value thereafter in the objective function.