👔 Hierarchical Optimization with Gurobi

One of the first technology choices to make when setting up an optimization stack is which modeling interface to use. Even if we restrict our choices to Python interfaces for MIP modeling, there are lots of options to consider.

If you use a specific solver, you can opt for its native Python interface. Examples include libraries like gurobipy, Fusion, highspy, or PySCIPOpt. This approach provides access to important solver-specific features such as lazy constraints, heuristics, and various solver settings. However, it can also lock you into a solver before ready for that.

You can also choose a modeling API that targets multiple solvers. In the Python ecosystem. These are libraries like amplpy, Pyomo, PyOptInterface, and linopy. These interfaces target multiple solver backends (both open source and commercial) and provide a subset of the functionality of each. Since they make it easy to switch between solvers, this is usually where I start.¹

Hierarchical assignment

However, there are plenty of times when solver-specific APIs are useful, or even critical. One example is hierarchical optimization. This is a simple technique for managing trade-offs between multiple objectives in a problem. Let’s look at an example.

Imagine we are assigning in-home health care workers ($w \in W$) to patients ($p \in P$). For simplicity, let’s say we have $n$ workers and $n$ patients, and we are assigning them one-to-one. Each worker has a given cost ($c_{wp}$) of assignment to each patient, which may reflect something like the travel time to get to them. We want to assign each worker to exactly one patient while minimizing the overall cost.

Model

So far, what we have is a simple linear sum assignment problem.

$$ \begin{align*} & \text{min} && z = \sum_{wp} c_{wp} x_{wp} \\ & \text{s.t.} && \sum_w x_{wp} = 1 && \forall \quad p \in P \\ & && \sum_p x_{wp} = 1 && \forall \quad w \in W \\ & && x \in \{0,1\}^{|W \times P|} \end{align*} $$

Solving this model gives us the minimum cost assignment. That’s all well and good, but now say we have a secondary objective of maximizing affinity of workers to patients ($a_{wp}$). That is, we want to prefer assignments that increase overall affinity while still minimizing cost. This is actually a common goal in health care scheduling: if possible, send the same worker to a given patient that you usually send.

Hierarchical optimization gives us a simple way to solve this problem. First, we optimize the model as stated above. This gives us an optimal objective value $z^*$. Then we re-solve the same optimization model, while constraining the cost to be $z^*$ and using the secondary objective function. This says to the optimizer, “improve the affinity as much as you can, but keep the cost optimal.”

$$ \begin{align*} & \text{max} && w = \sum_{wp} a_{wp} x_{wp} \\ & \text{s.t.} && \sum_{wp} c_{wp} x_{wp} \le z^* \\ & && \sum_w x_{wp} = 1 && \forall \quad p \in P \\ & && \sum_p x_{wp} = 1 && \forall \quad w \in W \\ & && x \in \{0,1\}^{|W \times P|} \end{align*} $$

From here, the natural question becomes: what if we trade off some cost for affinity? If we’re willing to increase cost by some percentage, how much more affinity do we get? We can do this by setting a constant $\alpha \ge 0$ and solving the model a number of times.²

$$ \begin{align*} & \text{max} && w = \sum_{wp} a_{wp} x_{wp} \\ & \text{s.t.} && \sum_{wp} c_{wp} x_{wp} \le (1 + \alpha) z^* \\ & && \sum_w x_{wp} = 1 && \forall \quad p \in P \\ & && \sum_p x_{wp} = 1 && \forall \quad w \in W \\ & && x \in \{0,1\}^{|W \times P|} \end{align*} $$

For example, if $\alpha = 0.05$, then we’re willing to accept a 5% increase in overall cost to improve affinity. Setting different values of $\alpha$ lets us explore the space of that trade-off and its impact on cost and affinity.

Once we solve this and get the optimal affinity ($w^*$), we should re-optimize for the primary objective again while constraining the secondary one.

$$ \begin{align*} & \text{min} && \sum_{wp} c_{wp} x_{wp} \\ & \text{s.t.} && \sum_{wp} a_{wp} x_{wp} \ge w^* \\ & && \sum_w x_{wp} = 1 && \forall \quad p \in P \\ & && \sum_p x_{wp} = 1 && \forall \quad w \in W \\ & && x \in \{0,1\}^{|W \times P|} \end{align*} $$

Code

So the math looks reasonable. How do we implement it? If we have a Gurobi license, we can use its built-in facilities for multiobjective optimization. This means that, instead solving a model multiple times and adding constraints to keep cost within $\alpha$ of its optimal value, we can create a single model that does all of this for us.

Assume we have input data which looks like this.

{
    "cost": [
        [10, 20, ...],
        [30, 40, ...],
        ...
    ],
    "affinity": [
        [25, 15, ...],
        [35, 25, ...],
        ...
    ]
}

We start with a simple assignment problem formulation.

import gurobipy as gp

n = len(data["cost"])
workers = range(n)
patients = range(n)

m = gp.Model()
m.ModelSense = gp.GRB.MINIMIZE

# x[w,p] = 1 if worker w is assigned to patient p.
x = m.addVars(n, n, vtype=gp.GRB.BINARY)

for i in range(n):
    # Each worker is assigned to one patient.
    m.addConstr(gp.quicksum(x[i, p] for p in patients) == 1)

    # Each patient is assigned one worker.
    m.addConstr(gp.quicksum(x[w, i] for w in workers) == 1)

We add primary and secondary objectives, and call optimize. The objectives are solved in descending order of the priority flag for Model.setObjectiveN. reltol allows us to degrade the primary objective by some amount (e.g. 5%) to improve the secondary objective.

One catch is that the model only has one objective sense. Since we are minimizing the primary objective, we give the secondary objective a weight of -1 in order to maximize it.

from itertools import product

# Primary objective: minimize cost.
z = (data["cost"][w][p] * x[w, p] for w, p in product(workers, patients))
m.setObjectiveN(expr=gp.quicksum(z), index=0, name="cost", priority=1, reltol=alpha)

# Secondary objective: maximize affinity. Since the model sense is minimize,
# we negate the secondary objective in order to maximize it.
w = (data["affinity"][w][p] * x[w, p] for w, p in product(workers, patients))
m.setObjectiveN(
    expr=gp.quicksum(w), index=1, name="affinity", priority=0, weight=-1
)

m.optimize()

Then we use this magic syntax to pull out the optimal cost and affinity.

m.params.ObjNumber = 0
cost = m.ObjNVal

m.params.ObjNumber = 1
affinity = m.ObjNVal

Results

If we solve this in a loop with alpha values from 0 to 1 in increments of 0.05, we can plot the trade-off between cost and affinity. Going from $\alpha = 0$ to $\alpha = 0.05$ or $\alpha = 0.1$ gives a pretty sizable improvement in affinity. After that, the return starts to gradually level off. This allows us to make a more informed choice about these two objectives.

Resources

• generate.py generates input data
• input-100x100.json contains input data
• model.py hierarchical objectives Gurobi model