I'm solving this optimization problem where I need to figure out how many distribution centers I need to open in order to fulfill the demand of the 12 company facilities, while minimizing the transportation costs. The transportation costs are simply the distance between the distribution centers times the cost per mile, however in this problem, the cost per mile is one dollar. I have 5 choices which are Boston, Nashua, Providence, Springfield and Worcester, these 5 are part of the 12 company facilities.
I solved the problem and got the correct answer but then I tried to add two constraints to the same code and the answer I get is incorrect. The two other constraints are that the average distance from the distribution centers (DC) to the other facilities (customer) has to be less than 60 miles; and the second constraint is that the percentage of customers within 50 miles has to be greater than 80% (0.8). I know the answer to this problem, the cost has to be $66,781 dollars, the average customer distance is 15 miles and the percentage of customers within 50 miles is 90%. The output of my code is that the cost is $66289 dollars, the average customer distance is 15.36 miles and the percentage of customers within 50 miles is 179%, which doesn't make sense.
Can you help me figure out why am I getting an odd output? Thanks in advance.
from gekko import GEKKO
import numpy as np
import pandas as pd
import math
m = GEKKO(remote=False) #So that it solves the problem locally
m.options.SOLVER = 1 #MILP
varx = [[0 for col in range(12)] for row in range(5)] #Creates an empty list
for i in range (5):
for j in range (12):
varx[i][j] = m.Var(lb = 0, integer = True)
varx = np.array(varx)
varxt = np.transpose(varx)
vary = np.empty([]) #Creates an empty array
for i in range(5):
vary = np.append(vary, m.Var(lb = 0, ub = 1, integer = True)) #Yes/No variables
vary = vary[1:13]
dists = np.array([[0 , 93, 69, 98, 55, 37, 128, 95, 62, 42, 82, 34], #Boston
[37, 65, 33, 103, 20, 0, 137, 113, 48, 72, 79, 41], #Nashua
[42, 106, 105, 73, 92, 72, 94, 57, 104, 0, 68, 38], #Providence
[82, 59, 101, 27, 93, 79, 63, 57, 127, 68, 0, 47], #Springfield
[34, 68, 72, 66, 60, 41, 98, 71, 85, 38, 47, 0]]) #Worcester
max_dist = 60 #Max average distance (in miles)
min_pct = 0.8 #Min percent of demand within 50 miles
aij = np.zeros((5, 12)) #Creates an empty array
for i in range (5):
for j in range (12):
if dists[i][j] <= 50:
aij[i][j] = 1
else:
aij[i][j] = 0 #Creates a 0s and 1s array. If the distance to a costumer
#is less than 50, then the matrix element is 1, it is zero
#otherwise
dem_consts = np.array([425, 12, 43, 125, 110, 86, 129, 28, 66, 320, 220, 182])
fixd_cost = 10000
sum1 = np.sum(np.multiply(varx, dists))
sum2 = np.sum(vary*fixd_cost)
z = sum1 + sum2
tot_dem = np.sum(dem_consts)
M = tot_dem
m.Minimize(z)
for i in range(12):
m.Equation(np.sum(varxt[i, :]) >= dem_consts[i]) #Demand constraints
for i in range(5):
m.Equation(np.sum(varx[i, :]) <= 2000) #Capacity constraints
m.Equation(np.sum(varx[i, :]) <= M*vary[i]) #Enforces 0 or 1 value
m.Equation(np.sum(vary[:]) >= 1)
di_sum = np.sum(np.multiply(varx, dists))
di_sumw = di_sum/ tot_dem
m.Equation(di_sumw <= max_dist) #Average (demand) weighted distance from DC to customer
a_sum = np.sum(np.multiply(varx, aij))
a_sumw = a_sum/tot_dem
m.Equation(a_sumw >= min_pct) #Percent of demand that is within 50 miles
m.solve(disp = False)
p1 = np.zeros((5, 12))
for i in range (5):
for j in range (12):
p1[i][j] = varx[i][j].value[0]
p1t = np.transpose(p1)
p2 = np.zeros((5, ))
for i in range(5):
p2[i] = vary[i].value[0]
mad1 = np.sum(np.multiply(p1, dists))
mad2 = mad1/tot_dem
mpi1 = np.sum(np.multiply(p1, aij))
mpi2 = mpi1/tot_dem
tot1 = np.sum(np.multiply(p1, dists))
tot2 = np.sum(p2)*fixd_cost
tot = tot1 + tot2
print('The minimum cost is:' +str(tot))
print('Average customer distance:' +str(mad2))
print('Percent of customers <= 50 miles:' +str(mpi2))
dc = np.array(['Boston', 'Nashua', 'Providence', 'Springfield', 'Worcester'])
cities = ['Boston', 'Brattleboro', 'Concord', 'Hartford', 'Manchester', 'Nashua',
'New Haven', 'New London', 'Portsmouth', 'Providence', 'Springfield', 'Worcester']
data = {cities[0]: p1t[0], cities[1]: p1t[1], cities[2]: p1t[2], cities[3]: p1t[3],
cities[4]: p1t[4], cities[5]: p1t[5], cities[6]: p1t[6], cities[7]: p1t[7],
cities[8]: p1t[8], cities[9]: p1t[9], cities[10]: p1t[10], cities[11]: p1t[11]}
df = pd.DataFrame(data, index = dc)
df
