How to dynamically compute nurse walking distance based on assignments in a nurse-patient scheduler model in Pyomo?

I am writing a Nurse-Patient assignment algorithm in Pyomo, but I am having trouble dynamically computing the walking distance between assigned patient beds for each nurse. The relevant sets and variables are below with a small example. The variables are filled in with the values that I would want at the end of solving.

#Sets
model.PatientIDs  = {0, 1, 2}
model.NurseIDs = {a, b}

#Indexed parameter (indexed by (patient1, patient2))
# Stores the distance between each patient bed  
model.Bed_Distances = {(0, 0) : 0, (0, 1) : 6.4, (0, 2) : 7.2, (1, 0) : 6.4, (1, 1) : 0, (1, 2) : 1.9, (2, 1) : 1.9), (2, 0) : 7.2}

# Indexed variable (indexed by (patient, nurse))
model.Assignments = {(0, a): 1, (0, b): 0, (1, a) : 0, (1, b) : 1, (2, a) : 1, (2, b) : 0}

# Indexed variable (indexed by nurse)
# Keeps track of the distance between assigned patient beds for each nurse
model.Nurse_Distances = {a : 7.2, b : 0}

I'm having trouble computing the nurse distances dynamically as the model is being solved (since it is dependent on the patients assigned to each nurse). Since the model.ASSIGNMENTS decision variable is represented with binary 0 or 1, I've been using this rule to compute the nurse distances:

def nurse_distance(self, nurse):
    return model.NURSE_DISTANCE[nurse] == sum([model.ASSIGNMENTS[p1, nurse] * model.ASSIGNMENTS[p2, nurse] * model.BED_DISTANCES[p1, p2] for p1 in model.PATIENTS for p2 in model.PATIENTS])

model.DISTANCE_CONSTRAINT = pe.Constraint(model.NURSES, rule = nurse_distance)

After adding this constraint, my model runs forever, or when I run it with only 2 or 3 patients it finishes relatively quickly and gives me an error that looks like this:

ERROR: Error parsing NEOS solution file NEOS log: Job 11992161 dispatched password: iGYfJtMj ---------- Begin Solver Output ----------packages/pyomo/opt/plugins/sol.py", line 87, in _load raise ValueError("no Options line found") ValueError: no Options line found

Is there a better way to keep track of a variable that is dependent on the value of another variable in Pyomo (nurse bed distances is dependent on the assignments made between nurses and patients)? If not, is there a more computationally efficient way to compute it within the constraint I've written without using "if then" statements?

Solution

You are creating a non-linear model by multiplying your decision variables together, which is probably not desired...

Because you have an "and" condition here (that the same nurse is assigned to p1 and p2) you almost certainly need another binary "helper" variable to track whether both are true, then use that. Of course, you will need a linking constraint between the individual assignments and the paired assignment. See the example below... Note that "Todd" gets off easy here because he only has 1 patient.

Code:

# Nurse Patrol

import pyomo.environ as pyo

# data
PatientIDs  = [0, 1, 2]
NurseIDs = ['Todd', 'Cindy']
Bed_Distances = {(0, 0) : 0, (0, 1) : 6.4, (0, 2) : 7.2, (1, 0) : 6.4, (1, 1) : 0, (1, 2) : 1.9, (2, 1) : 1.9, (2, 0) : 7.2}
patient_pairings = [(p1, p2) for p1 in PatientIDs for p2 in PatientIDs if p1 != p2]
# model
m = pyo.ConcreteModel('Nurse Patrol')
# put the stuff above into model parameters.  You can use them "raw" but this will help you T/S the model

# SETS
m.N = pyo.Set(initialize=NurseIDs)
m.P = pyo.Set(initialize=PatientIDs)
m.PP = pyo.Set(initialize=patient_pairings)

# PARAMS
m.distance = pyo.Param(m.P, m.P, initialize=Bed_Distances)

# VARIABLES
m.assign =          pyo.Var(m.N, m.P, domain=pyo.Binary)
m.pair_assigned =   pyo.Var(m.N, m.PP, domain=pyo.Binary)   # true if BOTH are assigned to nurse

# OBJ:  Minimize the total of all nurses' travels
m.obj = pyo.Objective(expr=sum(m.pair_assigned[n, p1, p2]*m.distance[p1, p2] 
                            for (p1, p2) in m.PP
                            for n in m.N))
# CONSTRAINTS
# Cover all patients
def covered(m, patient):
    return sum(m.assign[n, patient] for n in m.N) >= 1
m.C1 = pyo.Constraint(m.P, rule=covered)

# link the assignment to the paired assignment
def paired(m, n, p1, p2):
    return m.pair_assigned[n, p1, p2] >= m.assign[n, p1] + m.assign[n, p2] - 1
m.C2 = pyo.Constraint(m.N, m.PP, rule=paired)

result = pyo.SolverFactory('glpk').solve(m)
print(result)
m.assign.display()

for n in m.N:
    tot_dist = sum(m.pair_assigned[n, p1, p2]*m.distance[p1, p2] for (p1, p2) in m.PP)
    print(f'nurse {n} distance is: {pyo.value(tot_dist)}')

Yields:

Problem: 
- Name: unknown
  Lower bound: 3.8
  Upper bound: 3.8
  Number of objectives: 1
  Number of constraints: 16
  Number of variables: 19
  Number of nonzeros: 43
  Sense: minimize
Solver: 
- Status: ok
  Termination condition: optimal
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 3
      Number of created subproblems: 3
  Error rc: 0
  Time: 0.007218837738037109
Solution: 
- number of solutions: 0
  number of solutions displayed: 0

assign : Size=6, Index=assign_index
    Key          : Lower : Value : Upper : Fixed : Stale : Domain
    ('Cindy', 0) :     0 :   0.0 :     1 : False : False : Binary
    ('Cindy', 1) :     0 :   1.0 :     1 : False : False : Binary
    ('Cindy', 2) :     0 :   1.0 :     1 : False : False : Binary
     ('Todd', 0) :     0 :   1.0 :     1 : False : False : Binary
     ('Todd', 1) :     0 :   0.0 :     1 : False : False : Binary
     ('Todd', 2) :     0 :   0.0 :     1 : False : False : Binary
nurse Todd distance is: 0.0
nurse Cindy distance is: 3.8
[Finished in 558ms]