Implementing rule-based constraints featuring variables on different levels

[FranDeLio]

Hello Pao devs! First of all thank you for your awesome work!

I wanted to ask a question about how to correctly build the following bilevel optimization program. More specifically, I am struggling with the first constraint, as it features variables from both the model and the submodel.

$$\begin{align} %r0 & \underset{z_{j}}{\text{maximize}}~\underset{y_{ij}}{\text{minimize}} ~~ \sum_{i \in \mathcal{I}} \sum_{j \in \mathcal{J}} c_{ij}y_{ij}\\ &\\ s.t. \quad & \sum_{i \in \mathcal{I}} y_{ij} \leq z_j~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\forall j \in \mathcal{J}\\ & \sum_{j \in \mathcal{J}} y_{ij} \leq 1 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \forall i \in \mathcal{I}\\ %r5 & \sum_{i \in \mathcal{I}} \sum_{j \in \mathcal{J}} y_{cp} = min(|I|,|J|)\\ & y_{ij} \in \{0,1\} ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \forall (i,j) \in \mathcal{I} \times \mathcal{J}\\ & z_{j} \in \{0,1\} ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \forall (i,j) \in \mathcal{J}\\ \end{align}$$

If I initialize the constraint within the submodel e.g

model.submodel.each_taxi_serves_one_user_max = pe.Constraint(
            model.users, rule=each_taxi_serves_one_user_max
        )

File "/home/franciscodelio/francisco_python_projects/Dynamic_Assignment_Problem/static_optimization/solvers.py", line 311, in find_optimal_solution
    solver.solve(self.model, mip_solver="cbc")
  File "/home/franciscodelio/.cache/pypoetry/virtualenvs/dynamic-assignment-problem-ASJTY80e-py3.10/lib/python3.10/site-packages/pao/pyomo/solvers/mpr_solvers.py", line 36, in solve
    return super().solve(model, **options)
  File "/home/franciscodelio/.cache/pypoetry/virtualenvs/dynamic-assignment-problem-ASJTY80e-py3.10/lib/python3.10/site-packages/pao/pyomo/solver.py", line 55, in solve
    mp, soln_manager = convert_pyomo2MultilevelProblem(model, inequalities=True)
  File "/home/franciscodelio/.cache/pypoetry/virtualenvs/dynamic-assignment-problem-ASJTY80e-py3.10/lib/python3.10/site-packages/pao/pyomo/convert.py", line 554, in convert_pyomo2MultilevelProblem
    node.initialize_level(L, inequalities, var, vidmap, levelmap)
  File "/home/franciscodelio/.cache/pypoetry/virtualenvs/dynamic-assignment-problem-ASJTY80e-py3.10/lib/python3.10/site-packages/pao/pyomo/convert.py", line 204, in initialize_level
    t, nid, j = vidmap[vid]
KeyError: 140037957243440

then I can not access the variable in the model $z$ and I get the error above, and if I initialize the constraint from the model e.g

model.each_taxi_serves_one_user_max = pe.Constraint(
            model.users, rule=each_taxi_serves_one_user_max
        )

the program executes with no apparent errors, however the solutions provided violate the 1st constraint.

Hope you can help! Leaving a reproducible example below where the first constraint is initialized from the model (case 2).

from scipy.spatial import distance
import numpy as np
import pandas as pd
import pyomo.environ as pe
import pyomo.opt as po
import pao


"""Initializing dummy variables and useful functions"""
def generate_2d_coordinate():
    return tuple(np.random.uniform(low=0.0, high=100.0, size=2))

def all_unique(lst):
  return len(lst) == len(set(lst))

available_taxis_ids = [0,1,2,3,4]
users_ids = [0,1,2,3,4,6,7]

n_assignments = min(len(available_taxis_ids),len(users_ids))

taxi_coordinates = [
    generate_2d_coordinate() for _ in available_taxis_ids
]
destination_coordinates = [
    generate_2d_coordinate() for _ in users_ids
]
cost_matrix = distance.cdist(
    taxi_coordinates, destination_coordinates, "euclidean"
)

assignment_cost = (
    pd.DataFrame(
        cost_matrix, index=available_taxis_ids, columns=users_ids
    )
    .stack()
    .to_dict()
)

"""Build and solves the Pyomo model for the MILP problem."""
model = pe.ConcreteModel()

model.taxis = pe.Set(initialize=available_taxis_ids)
model.users = pe.Set(initialize=users_ids)
model.z = pe.Var(model.users, domain=pe.NonNegativeReals)

model.submodel = pao.pyomo.SubModel(fixed=[model.z])
model.submodel.taxis = pe.Set(initialize=available_taxis_ids)
model.submodel.users = pe.Set(initialize=users_ids)

model.submodel.taxis = pe.Set(initialize=available_taxis_ids)
model.submodel.users = pe.Set(initialize=users_ids)
model.submodel.assignment_cost = pe.Param(
    model.submodel.taxis, model.submodel.users, initialize=assignment_cost
)
model.submodel.y = pe.Var(model.submodel.taxis, model.submodel.users, domain=pe.NonNegativeReals)

expression = sum(
    model.submodel.assignment_cost[c, k] * model.submodel.y[c, k]
    for c in model.submodel.taxis
    for k in model.submodel.users
)

model.obj = pe.Objective(sense=pe.maximize, expr=expression)
model.submodel.obj = pe.Objective(sense=pe.minimize, expr=expression)

def each_taxi_serves_one_user_max(model, k):
    # assign at most one taxi c to every user.
    constraint = sum(model.submodel.y[c, k] for c in model.submodel.taxis) <= model.z[k]
    return constraint

def each_user_is_served_by_one_taxi_max(submodel, c):
    # a taxi c can only be assign to one given user k.
    constraint = sum(submodel.y[c, k] for k in submodel.users) <= 1
    return constraint

def maximal_injective_assignment(submodel):
    # a minimum of min(|Taxis|,|Users|) matches ought to be made.
    constraint = (
        sum(submodel.y[c, k] for k in submodel.users for c in submodel.taxis)
        == n_assignments
    )
    return constraint

model.each_taxi_serves_one_user_max = pe.Constraint(
    model.users, rule=each_taxi_serves_one_user_max
)

model.submodel.each_user_is_served_by_one_taxi_max = pe.Constraint(
    model.submodel.taxis, rule=each_user_is_served_by_one_taxi_max
)
model.submodel.maximal_injective_assignment = pe.Constraint(
    rule=maximal_injective_assignment
)

with pao.pyomo.Solver('pao.pyomo.FA') as solver:
    solver.solve(model, mip_solver="cbc")

"""Processing Results"""

output_df = (
            pd.Series(model.submodel.y.extract_values())
            .reset_index()
            .rename(columns={"level_0": "taxis", "level_1": "users", 0: "y"})
        )

solution_is_binary = set(output_df.y.unique()).issubset({0, 1})
assert solution_is_binary

solution_df = output_df.loc[lambda x: x.y == 1, ["users", "taxis"]].reset_index(
            drop=True
        )

print(solution_df)
assert all_unique(solution_df.users)

users  taxis
0      0      0
1      2      1
2      6      2
3      0      3
4      7      4

AssertionError                            

[FranDeLio]

Update: I was able to solve this problem using Julia's Bilevel.JuMP.