[Python] How to solve this linear programming problem?

Consider a set of city districts . There is a candidate for building a medical emergency service station in every city district. The cost of building an ambulance station is given by the parameter and the total budget for building stations throughout the city is 50 million. Furthermore, if a station is built in city district 2, it must not be station built in city district 6. Determine in which city districts to build ambulance stations so that the maximum travel time is minimized. What is the optimal maximum travel time? The cost of building stations in the city district ∈ , are: [32,20,25,30,40,29] The city districts and the driving distances between them are interpreted in the following graph:

([[ 0, 4, 12, 27, 25, 58],
[ 4, 0, 24, 16, 29, 38],
[12, 24, 0, 31, 14, 30],
[27, 16, 31, 0, 21, 8],
[25, 29, 14, 21, 0, 11],
[58, 38, 30, 8, 11, 0]])


The right answer is to build stations in districts 2 and 4, and optimal travel time is 24.

I tried to solve it, and here is my code, but the answer is not right. It has to be completet with PULP

from pulp import *
import numpy as np

# the cost of building of the stations
cost1 = [32,20,25,30,40,29]

# the city districts and the driving distances between them
d = np.array([[0,4,12,27,25,58],[4,0,24,16,29,38],[12,24,0,31,14,30],[27,16,31,0,21,8], [25,29,14,21,0,11],[58,38,30,8,11,0]])

# set of city districts
cities1 = [1, 2, 3, 4, 5, 6]
cities2 = [1, 2, 3, 4, 5, 6]

# costs
n = dict(zip(cities1,cost1))

# dict of every distance between the cities
a1 = dict(zip(cities1, zip(*d)))

a = {key: dict(enumerate(value, 1)) for key, value in a1.items()}


# model
model = LpProblem("Set_covering", LpMinimize)

# shortcut introduction
shortcut = [(i,j) for i in cities1 for j in cities2]

# declaration of variables
x = LpVariable.dicts("x",(cities1),0, cat = 'Binary')
maximum = LpVariable("Maximum", 0, cat = 'Continuous')

# purpose function
model += maximum

# restrictive conditions
for (i,j) in shortcut:
model += maximum >= x*a[j]
for j in cities1:
model += lpSum(x*a[j] for i in cities2) >= 1

model += lpSum(n*x for i in cities1) <=50
model += lpSum(x[2]+x[6]) <= 1


# List of results
model.solve()
print("Status:", LpStatus[model.status])

for v in model.variables():
if v.varValue > 0:
print(v.name, "=", v.varValue)

print("The total travel time to all places from both stations is: ", value(model.objective))`

Continue reading...