infeasible results without t0 constraints, errors with using t 0 constraints

Hi,

I wonder if someone could guide me on where I am going wrong here.

If I remove the constraint for each time period, the problem is solved with infeasible, unbounded.

Once I added the time constraint, I get the following error:

import itertools
from itertools import product
from math import sqrt
import pandas as pd
import gurobipy as gp
from gurobipy import GRB
from gurobipy import quicksum


m = gp.Model('VCStrategy')

#minimum amount of dose made by one reactor per day (one minute run per day)
minDosePerR = 600

#max number of hours every reactor can run for (8 hours a day)
maxRunTime = 8*60

#county population
CountyPopulation = 7000

centers = [(0,0), (0,1), (0,2), (1,0), (1,1), (1,2), (2,0), (2,1), (2,2)] #this needs to reflect the db
facilities = [(0,0), (0,1), (0,2)] #this needs to reflect the db

# Indices for the facilities
facilitiesIndex = [*range(1,3)] #this needs to reflect the count from db

# Indices for the counties / later to be post codes/ as in GPs
countiesIndex = [*range(1,9)] #this needs to reflect count the db

#time period - 14 days
period = [*range(0,13)] 

#nurse number of admissions per day
nurse_ratePerDay = 28 

#max number of nurses in each clinic to administrate vaccines
nurse_max = 8 

#cost to set up one reactor + storage staff and the rest for 14 days
costPerFacilitySetup = [10000,10000,10000] #this needs to reflect the db for count of hospitals provided 10,000

#cost of travel per mile by vehicle
costPerMile = 10

#cost per day per nurse
costPerNurse = 105

capacityTOStorePerGP = 399
capacityTOStorePerHospital = 2000

storageCostPerDoseGP = 0.90
storageCostPerDoseHospital = 0.79


budget = 200000

# Cartesian product of counties and facilities
fc = [] 
for f in facilitiesIndex:
    for c in countiesIndex:
        tp = f,c
        fc.append(tp)

# Cartesian product of facilites and centers and time period
fct = []
for f in facilitiesIndex:
    for c in countiesIndex:
        for t in period:
            tp = f,c,t
            fct.append(tp)

# Cartesian product of facilities and time
ft = []
for f in facilitiesIndex:
    for t in period:
        tp = f,t
        ft.append(ft)

# Cartesian product of centers and time
ct = []

for c in countiesIndex:
    for t in period:
        tp = c,t
        ct.append(tp)
        
        
# This function determines the Euclidean distance between a facility and customer sites.
def compute_distance(loc1, loc2):
    dx = loc1[0] - loc2[0]
    dy = loc1[1] - loc2[1]
    return sqrt(dx*dx + dy*dy)

# Compute key parameters of MIP model formulation

num_facilities = len(facilities)
num_customers = len(centers)
cartesian_prod_fc = list(product(range(num_facilities), range(num_customers)))
cartesian_prod_fct = list(product(range(num_facilities), range(num_customers), (period)))
cartesian_prod_ct = list(product(range(num_customers), (period)))
cartesian_prod_ft = list(product(range(num_facilities), (period)))


select = m.addVars(num_facilities, vtype=GRB.BINARY, name='Select')

# Compute shipping costs

shipping_cost = {(f,c): costPerMile*compute_distance(facilities[f], centers[c]) for f, c in cartesian_prod_fc}

#number of reactors at facility f at period t
rn = m.addVars(cartesian_prod_ft, vtype=GRB.CONTINUOUS, lb = 0,  name="reactor_{}".format((f,t) for f,t in cartesian_prod_ft))

#number of doses at facility f at period t
d = m.addVars(cartesian_prod_ft, vtype=GRB.CONTINUOUS, lb = 0,  name="doses_{}".format((f,t) for f,t in cartesian_prod_ft))

#DurationOfRunPerReacorInTimePeriod - #minimum running time is 1 minute lower bound
rnt = m.addVars(cartesian_prod_ft, vtype = GRB.CONTINUOUS, lb = 1, name = "reactorRunTime_{}".format((f,t) for f,t in cartesian_prod_ft))

#rn*rnt*600 = total doses per day

#number of reactors * duration of run * 600

#number of allocated doses from facility f for center c at period t
alloc = m.addVars(cartesian_prod_fct, vtype=GRB.CONTINUOUS, name="allocated_{}".format((f,c,t) for f,c,t in cartesian_prod_fct))

#shipped doses
s = m.addVars(cartesian_prod_fct, vtype=GRB.CONTINUOUS, name="Shipped_{}".format((f,c,t) for f,c,t in cartesian_prod_fct))

#administrated doses
adm = m.addVars(cartesian_prod_ct, vtype=GRB.CONTINUOUS, name="administrated_{}".format((c,t) for c,t in cartesian_prod_ct))

#doses stored at facility f in period 
fstored = m.addVars(cartesian_prod_ft, vtype=GRB.CONTINUOUS, name="fstored{}".format((f,t) for f,t in cartesian_prod_ft))

#doses stored at center c in period
cstored = m.addVars(cartesian_prod_ct, vtype=GRB.CONTINUOUS, name="cstored_{}".format((c,t) for c,t in cartesian_prod_ct))

#number of nurses at location c
nc =  m.addVars(cartesian_prod_ct, vtype = GRB.CONTINUOUS, lb = 0, name = "numberofNurses_{}".format((c,t) for c,t in cartesian_prod_ct))

#number of deliveries from factory to center in period
nd =  m.addVars(cartesian_prod_fct, vtype = GRB.CONTINUOUS, lb = 0, name = "deliveries_{}".format((f,c,t) for f,c,t in cartesian_prod_fct))                              




#List of Contraints

# m.addConstr(quicksum(d[f, t]) for f,t in cartesian_prod_ft == CountyPopulation)




#total of doses made in 14 days need to equal demand in county
m.addConstrs((d[f, t] == CountyPopulation for f,t in cartesian_prod_ft
                            if f != t), name='demandConstraints')

#max hour of operation per reactor per day is 480 minutes
m.addConstrs((rnt[f, t] <= maxRunTime for f,t in cartesian_prod_ft
                            if f != t), name='maxRunTimeConstraints')

#max number of nurses shouldn't be more than 8
m.addConstrs((nc[c, t] <= maxRunTime for c,t in cartesian_prod_ct
                            if c != t), name='maxNurseConstraints')


#minimum produceddosesin each period is 600 doses if the factory is selected - as in the hospital is the point of production

#if t=0, none of the centers have any excess doses
#total allocated doses from facility f to center c in time period t shouldn't be greater than center's capacity for the day which is 399

#if t!=0, centers have some remaining doses 
#total allocated doses from facility f to center c in time period t shouldn't be greater than center's capacity - remaining doses

for t in period:
        if t == 0:
            # Counties demand constraints
            m.addConstrs((alloc[f,c,t] <= capacityTOStorePerGP for f,c,t in fct), name="OverTimeAllocation_{}".format(f,c,t))
            m.addConstrs((fstored[f,t] <= capacityTOStorePerHospital for f,t in cartesian_prod_ft), name="facilityStorageCapacity")
            m.addConstrs((cstored[c,t] <= capacityTOStorePerGP for c,t in ct), name="GPStorageCapacity")
    
        else:
            #allocated doses needs to be less than the capacity to store minus the ones already in storage
            m.addConstr((gp.quicksum(alloc[f, c, t]  <= (capacityTOStorePerGP - (alloc[f,c,t-1] - adm[c,t-1])) for f,c,t in fct))
                         , name = "OverTimeAllocation_{}".format(f,c,t))
            
            #difference between what was allocated and how much of that was actually used
            tadm = (alloc[f,c,t-1] for f,c,t in fct - adm[c,t-1] for c,t in ct)
            
            #difference between center's capacity (so how many are now left over in storage) and how much of that was actually used
            mius = (capacityTOStorePerGP - tadm for f,c,t in fct)

#             m.addConstrs((alloc[f,c,t] <= mius for f,c,t in fct), name="capacityTOStorePerGP")
#             m.addConstrs((fstored[f,t] <= capacityTOStorePerHospital for f,t in cartesian_prod_ft), name="facilityStorageCapacity")
#             m.addConstrs((cstored[c,t] <= capacityTOStorePerGP for c,t in ct), name="GPStorageCapacity")
    
            m.addConstrs((fstored[f,t] <= capacityTOStorePerHospital - fstored[f,t-1]  for f,t in ft),name="facilityStorageCapacity")
            m.addConstrs((cstored[c,t] <= capacityTOStorePerGP - cstored[c,t-1] for c,t in ct),name="GPStorageCapacity")
 


#total cost of setting the factories, shipping doses to all locations, administrating it and storing them in both factory and center shouldn't be greater than our budget for the total of 14 days

#calculate the total for setting up the factories, shipping the doses and administrating it
#total cost involves number of reactors time their fixed fee + numberofdeliveries * shippingcost + numberofstaff*theirdaily rate
#this needs to reflect all 14 period

#not sure how to do it for the whole 14 day period

m.update()


#min total excess in all periods = Total of produced in all periods - (total of all allocated per clinic - total of all dispensed)
# integrate variables and constraints
    

# for Objective function
totalPro = quicksum(d[f,t] for f,t in cartesian_prod_ft)
totalAllocated = quicksum(alloc[f,c,t] for f,c,t in fct)
totalAdm = quicksum(adm[c,t] for c,t in ct)


obj = (totalPro - (totalAllocated - totalAdm))
m.setObjective(obj,GRB.MINIMIZE)

# start optimization
m.optimize()

<ipython-input-34-2395b4ad79e3> in <genexpr>(.0)
    180         else:
    181             #allocated doses needs to be less than the capacity to store minus the ones already in storage
--> 182             m.addConstr((gp.quicksum(alloc[f, c, t]  <= (capacityTOStorePerGP - (alloc[f,c,t-1] - adm[c,t-1])) for f,c,t in fct))
    183                          , name = "OverTimeAllocation_{}".format(f,c,t))
    184 

KeyError: (1, 1, -1)