Infeasible model with VRPTW

model = Model("Drone Optimization")

N = 5

D = 2

q = {i: 1 for i in range(N)}

C = 10

flight_time_matrix = [

[0, 15, 20, 35, 50], # hospital 0

[15, 0, 10, 25, 40], # hospital 1

[20, 10, 0, 15, 30], # hospital 2

[35, 25, 15, 0, 20], # hospital 3

[50, 40, 30, 20, 0] # hospital 4

]

time_windows = {

0: (0, 1440), # Depot: always open

1: (120, 240), # H 1: 02:00 - 04:00

2: (300, 420), # H 2: 05:00 - 07:00

3: (480, 600), # H 3: 08:00 - 10:00

4: (660, 780), # H 4: 11:00 - 13:00

}

x = model.addVars(N, N, D, vtype=GRB.BINARY, name="x") #Drone k flies from i to j

y = model.addVars(N, D, vtype=GRB.BINARY, name="y") # Hospital j visted by drone k

arrival_time = model.addVars(D, N, vtype=GRB.CONTINUOUS, name="arrival_time")

# Minimize travel time

model.setObjective(quicksum(flight_time_matrix[i][j] * x[i, j, k]

for i in range(N) for j in range(N) for k in range(D) if i != j),

GRB.MINIMIZE)

# Constraints

# Every customer without a depot is visited once.

model.addConstrs((quicksum(y[j, k] for k in range(D)) == 1 for j in range(1,N)), name="visit_once")

# every drone visits a hospital if it was assigned

model.addConstrs(

(quicksum(x[i, j, k] for j in range(N) if j != i) == y[i, k])

foriinrange(1, N) forkinrange(D)

)

# every drone is doing a journey

model.addConstrs(

quicksum(x[i, j, k] foriinrange(N) forjinrange(N) ifi != j) >= 1

forkinrange(D)

)

# flow conservation

model.addConstrs(

quicksum(x[i, j, k] forjinrange(N) ifj != i) == quicksum(x[j, i, k] forjinrange(N) ifj != i)

forkinrange(D) foriinrange(N)

)

model.addConstrs(

quicksum(x[i, j, k] forjinrange(N) ifj != i) == y[i, k]

forkinrange(D) foriinrange(N)

)

# capacity constraint

model.addConstrs(

(quicksum(q[i] * y[i, k] for i in range(N)) <= C

forkinrange(D))

)

model.addConstrs(

quicksum(x[0, j, k] forjinrange(1, N)) == 1forkinrange(D) # Jede Drohne startet bei 0

)

model.addConstrs(

quicksum(x[i, 0, k] foriinrange(1, N)) == 1forkinrange(D) # Jede Drohne endet bei 0

)

#subtour elimination

u = model.addVars(range(N), vtype=GRB.INTEGER, lb=0, ub=N-1, name="u")

model.addConstrs(

u[i] - u[j] + (N-1) * x[i, j, k] <= N-2

forkinrange(D) foriinrange(1, N) forjinrange(1, N) ifi != j

)

model.addConstrs(u[0] == 0 for k in range(D))

M = 1000 # Große Zahl für die Modellierung der Bedingungen

# model.addConstrs(

# arrival_time[k, 0] == 0

# for k in range(D)

# )

# arival time for drone k flying from i to j

model.addConstrs(

arrival_time[k, j] >= arrival_time[k, i] + 60 + flight_time_matrix[i][j] - M * (1 - x[i, j, k])

forkinrange(D) # Für jede Drohne

foriinrange(N) # Für jedes Krankenhaus als Startpunkt

forjinrange(N) # Für jedes Krankenhaus als Zielpunkt

ifi != j# Nur für i != j, damit eine Drohne nicht von sich selbst wegfliegt

)

# z = model.addVars(D, N, N, vtype=GRB.CONTINUOUS, name="z") # Hilfsvariable für den maximalen Wert

# model.addConstrs(

# z[k, i, j] >= arrival_time[k, i] + 60 + flight_time_matrix[i][j] - M * (1 - x[i, j, k]) # Bedingung für den maximalen Wert

# for k in range(D)

# for i in range(N)

# for j in range(N)

# if i != j

# )

# # model.addConstrs(

# # arrival_time[k, j] >= z[k, i, j] # Ankunftszeit-Bedingung

# # for k in range(D)

# # for i in range(N)

# # for j in range(N)

# # if i != j

# # )

# # Optional: z[k, i, j] darf niemals negativ sein

# model.addConstrs(

# z[k, i, j] >= 0

# for k in range(D)

# for i in range(N)

# for j in range(N)

# if i != j

# )

# time window constraints

model.addConstrs(

(arrival_time[k, j] >= time_windows[j][0] for k in range(D) for j in range(1, N)), # earliest time

name="time_window_start"

)

model.addConstrs(

(arrival_time[k, j] <= time_windows[j][1] for k in range(D) for j in range(1, N)), # latest time

name="time_window_end"

)

# optimize

model.optimize()

# results

if model.status == GRB.OPTIMAL:

print("Optimale Lösung gefunden:")

forkinrange(D):

foriinrange(N):

forjinrange(N):

ifx[i, j, k].x > 0.5:

arrival_time_value = arrival_time[k, j].x

print(f"Fahrzeug {k}: {i} -> {j}, Ankunftszeit: {arrival_time_value}")

forkinrange(D):

forjinrange(N):

ify[j, k].x > 0.5:

print(f"Fahrzeug {k} besucht Kunde {j}")

else:

print("Keine optimale Lösung gefunden.")