The year is $1423$ and the storm is gathering around Lombaria! The ambitious Filippo Maria Visconti, Duke of Milan, have attacked you the fellow Venetians! You and your mercenary friends Tuna, Kaan, Ali, Kutup, O $$zgu $$n and Muhittin are stuck in Lombardy, Italy. You need to somehow save the people of Lombardy from the cruel Duke. You are currently located in Pavia, and you can rescue people to come with you along your journey in Lombardy. Due to the logistic reasons, each of you have the a capacity to rescue $82$ people during your route. Since the Duke is fast and furious, you only have time to visit $7$ cities in total.

But you are lucky, all of you are Medieval IE$/$OR specialists! So$,$ your job is to gather as many people as possible in order to collaborate and survive against the Milanese. In other words, you need to come up with a route that maximizes the number of people rescued along your route.

Rules

$$ From a city, you can only go to a neighboring city. $$ You need to start from Pavia.

$$ The route is finished once you visit the $7$th city.

The region of Lombardy has $12$ cities in total that can be seen in Figure $1.$

We also have some additional conditions. Ali knows someone named Sena in Milan, who must be rescued since she is the only medical scientist left in Lombardy. Moreover, Kutup knows a new companion of his, Bora, with a capacity to rescue $78$ more people, who is living right next to the Lake Como. O $$zgu $$n knows a person in Bergamo, Yi $$git$,$ who is among the population of Bergamo. O $$zgu $$n informs you that if Yi $$git is to come across with his $36$ friends in Cremona, Yi $$git is going to persuade them to leave your party and they will rush to the aid of the Venetians in the battlefield against the Milanese, so they will NOT be rescued in the end.

Solve this problem using Mixed Integer Linear Programming and report your solution The region of Lombardy has $12$ cities in total that can be seen in Figure $1.$

Figure $1$: Map of the region of Lombardy, Italy

Each city has a number of people residing in it as follows:after implementing it in gurobipy. Clearly indicate the visited cities in your reports.

For this question, I have the following code:

import gurobipy as gp

from gurobipy import GRB

# Cities & Populations

cities $=["$Pavia$",$"Lodi", "Milan", "Varese", "Como", "Monza$/$Brianza$",$ "Sondrio",

"Lecco","Cremona", "Mantua", "Brescia", "Bergamo"$]$

p $=[86,23,326,89,60,86,18,34,36,41,127,111]$

K $=[$

$[0,1,1,0,0,0,0,0,0,0,0,0],$

$[1,0,1,0,0,0,0,0,1,0,0,0],$

$[1,1,0,1,0,1,0,0,1,0,0,1],$

$[0,0,1,0,1,1,0,0,0,0,0,0],$

$[0,0,0,1,0,1,1,1,0,0,0,0],$

$[0,0,1,1,1,0,1,0,1,0,0,1],$

$[0,0,0,0,1,0,0,1,0,0,1,1],$

$[0,0,0,0,1,1,1,0,0,0,0,1],$

$[0,1,1,0,0,0,0,0,0,1,1,1],$

$[0,0,0,0,0,0,0,0,1,0,1,0],$

$[0,0,0,0,0,0,1,0,1,1,0,1],$

$[0,0,1,0,0,1,1,1,1,0,1,0]$

$]$

# Initialize model

m $=$ gp$.$Model$("$Italy$")$

# Decision variables

#X $=$ arc between cities i and j

#Y $=$ Visited city

#C $=$ Capacity of every person

#B $=$ Binary value that demonstrates whether bora is saved or not

x $=$ m$.$addVars$($len$($cities$),$ len$($cities$),$ vtype$=$GRB$.$BINARY, name$="$x$\_$ij$")$

y $=$ m$.$addVars$($len$($cities$),$ vtype$=$GRB$.$BINARY, name$="$y$\_$i$")$

c $=$ m$.$addVars$($range$(1,8),$ vtype$=$GRB$.$INTEGER, name$="$c$\_$i$")$

b $=$ m$.$addVar$($vtype$=$GRB$.$BINARY, name$=$"bora"$)$

ns $=$ m$.$addVar$($vtype$=$GRB$.$BINARY, name$=$"Not Saved"$)$

m$.$update$()$

# Objective Function

m$.$setObjective$($gp$.$quicksum$($c$[$i$]$ for i in range$(1,8))+78*$b $-$ ns$*$p$[8],$ GRB$.$MAXIMIZE$)$

# Each city is visited at most once

for j in range$($len$($cities$))$:

m$.$addConstr$($gp$.$quicksum$($K$[$i$][$j$]*$ x$[$i$,$ j$]$ for i in range$($len$($cities

# Constraints

#Starting from Pavia

m$.$addConstr$($y$[0]>=1)$

# Max. $7$ cities can be visited

m$.$addConstr$($gp$.$quicksum$($x$[$i$,$ j$]$ for i in range$($len$($cities$))$ for j in range$($len$($cities

for i in range $($len$($cities$))$:

for j in range $($len$($cities$))$:

m$.$addConstr$($K$[$i$][$j$]>=$ x$[$i$,$j$])$

#Model ends at most $7$th visit $(8$th city$)$

m$.$addConstr$($gp$.$quicksum$($y$[$i$]$ for i in range$($len$($cities

# Each arc can be used once

for i in range $($len$($cities$))$:

for j in range $($len$($cities$))$:

m$.$addConstr$($x$[$i$,$j$]+$x$[$j$,$i

# Capacity for each person

for j in range$(1,8)$:

m$.$addConstr$($c$[$j

# Sena, Bora & Yigit's special cases

m$.$addConstr$($y$[2]>=1)$

#m$.$addConstr$($y$[11]+$ y

if m$.$addConstr$($y$[11]==1)$:

ns $==1$

if m$.$addConstr$($y$[4]==1)$:

b $==1$

# To ensure that capacity does not exceeds population

m$.$addConstr$($gp$.$quicksum$($y$[$i$]*$ p$[$i$]$ for i in range$($len$($cities$)))>=$ gp$.$quicksum$($c$[$k$]$ for k in range$(1,8)))$

# To connect variables X & Y

for i in range $($len$($cities$))$:

for j in range $($len$($cities$))$:

if i$!=$j:

m$.$addConstr$($y$[$i$]+$ y$[$i$]>=2*$ x$[$i$,$j$])$

However, when the code is runned, the outcome is Pavia, Milan, Como, Bergamo. Since milan and como are not neighbor, it is impossible. Enhance the code according to that "neighbour" problem please.