import pandas as pd
import numpy as np
from scipy.spatial import distance
from mip.model import *
factories = pd.DataFrame ([
['F_01', 59.250276, 59.871183, 389],
['F_02', 84.739320, 14.179336, 409],
['F_03', 42.397937, 42.474530, 124],
['F_04', 19.539202, 13.714643, 70],
['F_05', 41.280669, 37.860993, 386],
['F_06', 37.159066, 41.353602, 196],
['F_07', 96.890453, 64.420010, 394],
['F_08', 86.267499, 81.662811, 365],
], columns=['fact', 'x', 'y', 'supply']).set_index ('fact')
shops = pd.DataFrame ([
['S_01', 13.490869, 73.269974, 200],
['S_02', 85.435435, 66.637250, 20],
['S_03', 28.578297, 8.997380, 320],
['S_04', 31.324145, 91.839907, 360],
['S_05', 40.338575, 15.487028, 360],
['S_06', 41.642451, 42.121572, 120],
['S_07', 53.983692, 20.950457, 360],
['S_08', 75.761895, 87.067552, 60],
['S_09', 81.836739, 36.799647, 80],
['S_10', 54.260517, 25.920108, 100],
['S_11', 67.918105, 68.108601, 340],
['S_12', 92.200710, 10.898110, 360],
['S_13', 19.966539, 39.046271, 60],
], columns=['shop', 'x', 'y', 'demand']).set_index ('shop')
unit_revenue = 1200
unit_production_cost = 200
unit_delivery_cost_per_km = 20
shop_rent_cost = 100000
N_facts = len (factories)
N_shops = len (shops)
to_df = lambda arr: pd.DataFrame (arr, index=shops.index, columns=factories.index)
factories_vs_shops = lambda fn: [[*(fn (f, s) for _, f in factories.iterrows ())] for _, s in shops.iterrows ()]
unit_delivery_cost = lambda f, s: (distance.euclidean ((f.x, f.y), (s.x, s.y)) * unit_delivery_cost_per_km)
unit_profit = lambda f, s: unit_revenue - (unit_production_cost + unit_delivery_cost (f, s))
unit_profits = factories_vs_shops (unit_profit)
max_supply = factories.supply.values
max_demand = shops.demand.values
max_total_supply = np.sum (max_supply)
def display_and_check_results (units, shops_open = np.ones (N_shops)):
supply = np.sum (units, axis=0)
demand = np.sum (units, axis=1)
df = to_df (units).replace (0.0, '')
shop_profits = to_df (units * unit_profits).sum (axis=1)
df['profit'] = shop_profits
df['rent'] = np.array (shops_open) * shop_rent_cost
df['total'] = shop_profits - df['rent']
df.loc[''] = None
df.loc['total'] = df.sum (axis=0)
display (df.fillna (''))
grand_total = df['total']['total']
print (f'Grand total: {grand_total}')
# All supply should be used
assert (max_total_supply == np.sum (supply))
# Supply & demand constraints should be satisfied!
assert (np.alltrue (np.round (supply) == max_supply) and np.alltrue (np.round (demand) <= max_demand))
# --------------------- SOLVE ----------------------------
m = Model (sense=MAXIMIZE, solver_name=CBC)
profits_coeffs = np.array ([[unit_profits[s][f] for f in range (N_facts)] for s in range (N_shops)]).flatten ()
# N×M matrix of decision variables (shop-wise)
supplies_s = [[*(m.add_var (var_type=CONTINUOUS, lb=0) for _ in range (N_facts))] for _ in range (N_shops)]
# The same but factory-wise
supplies_f = np.transpose (supplies_s)
# The same but all-flat
supplies = np.array (supplies_s).flatten ()
# Binary variables encoding open/closed state for shops
shops_open = [m.add_var (var_type=BINARY) for _ in range (N_shops)]
# Add linear constraints
for s in range (N_shops): m += xsum (supplies_s[s]) <= float (max_demand[s]) * shops_open[s]
for f in range (N_facts): m += xsum (supplies_f[f]) == float (max_supply[f])
# Objective function (what we are optimizing)
m.objective = (xsum (supplies[i] * profits_coeffs[i] for i in range (len (supplies))) -
xsum (shops_open[i] * shop_rent_cost for i in range (N_shops)))
# Run it!
display (m.optimize ())
display (m.objective_value)
display_and_check_results (np.array ([v.x for v in supplies]).reshape (N_shops, N_facts),
[v.x for v in shops_open])