python - PuLP: How to structure this factory coverage problem

Question

Struggling to formulate the objective function to the following MIP scenario in PuLP.

Each city is populated by N number of people.
Each factory location can facilitate a set of cities.
Minimize the number factories to open, such that the number of people facilitated is >= 4000

My main issue comes from the fact that different factories can service the same cities. So it isn't fair to sum the servicable population population of each factory and consider them seperately.

cities = ['London', 'Paris', 'Berlin', 'Amsterdam', 'Vienna', 'Prague']
factories = ['A', 'B', 'C', 'D']

city_populations = {'London': 898, 'Paris': 222, 'Berlin': 767, 'Amsterdam': 111, 'Vienna': 854, 'Prague': 908}

factories_service = {'A': ['London', 'Prague'], 'B': ['London', 'Paris', 'Vienna'], 'C': ['Amsterdam', 'Vienna', 'Prague'], 'D': ['London', 'Vienna', 'Prague']}

This is what I have at the moment but it is incorrect as it just picks the largest cities with no regard for population overlap.

prob = pl.LpProblem("Factory Coverage",pl.LpMinimize)
?
decision_vars = pl.LpVariable.dicts("Factories", factories, cat='Binary')
?
prob += pl.lpSum(decision_vars)
prob += pl.lpSum([sum([city_populations[x] for x in factories_service[i]])*decision_vars[i] for i in factories]) >= 4000
?
?
prob.solve()

Output:
Factories_A, 0.0
Factories_B, 1.0
Factories_C, 0.0
Factories_D, 1.0

Answer 1

The main problem is that you are double-counting the population. If two factories service the same city, and both are built, the population of that city should be accounted as serviced only once.

You can accomplish that by creating another set of binary variables for the cities, indicating if the city is serviced or not. And then add logic that determines that a city is serviced only if there is a factory that services it.

Another minor problem you have is that you don't have 4000 citizens in your entire model since sum(city_populations.values()) == 3760

This is an example of how you could solve this problem, though it is not the most efficient implementation:

prob = pl.LpProblem("Factory Coverage",pl.LpMinimize)

factories_vars = pl.LpVariable.dicts("Factories", factories, cat='Binary')
cities_vars = pl.LpVariable.dicts("Cities", cities, cat='Binary')

# Objective function
prob += pl.lpSum(factories_vars)
# Population served is greater than 4000
prob += pl.lpSum([city_populations[c]*cities_vars[c] for c in cities]) >= 2000
# City is served only if there is a factory that serves it
bigM = len(factories)
for c in cities:
    prob +=  bigM * cities_vars[c] >= pl.lpSum([factories_vars[f] for f in factories_vars if c in factories_service[f]])
    prob +=  cities_vars[c] <= pl.lpSum([factories_vars[f] for f in factories_vars if c in factories_service[f]])

prob.solve()

Result:

for f in factories_vars:
    print(f, factories_vars[f].value())
for c in cities_vars:
    print(c, cities_vars[c].value())

A 0.0
B 0.0
C 0.0
D 1.0
London 1.0
Paris 0.0
Berlin 0.0
Amsterdam 0.0
Vienna 1.0
Prague 1.0