设施位置 - 最小化为具有距离限制的客户提供服务的设施的算法

Question

例如，我有 1000 位客户位于不同纬度和经度的欧洲。我想找到可以为所有客户提供服务的最少数量的设施，受限于每个客户必须在 24 小时交付内得到服务（这里我使用从设施到客户的最大允许运输距离作为确保 24 小时交付的约束服务（距离是两个位置之间的直线，根据欧氏distance/straight线计算）。

因此，每个仓库只能为一定距离内的客户提供服务，例如600 公里，有什么算法可以帮助我找到为所有客户提供服务所需的最少设施数量，以及它们各自的经纬度。下面的附图中显示了一个示例。

查找最小仓库及其位置的示例

Answer 1

这属于设施位置问题的范畴。关于这些问题有相当丰富的文献。 p-center问题接近你想要的。

一些注意事项：

除了解决正式的数学优化模型外，通常还会使用启发式（和元启发式）。
这些距离是实际旅行时间的粗略近似值。这也意味着近似解可能已经足够好了。
除了找到为所有客户提供服务所需的最少数量的设施外，我们还可以通过最小化距离来优化位置。
纯“最小化设施数量”的数学规划模型可以表述为混合整数二次约束问题 (MIQCP)。这可以使用标准求解器（例如 Cplex 和 Gurobi）来解决。下面是我拼凑的一个例子：

通过 1000 个随机客户位置，我可以找到经过验证的最佳解决方案：

    ----     57 VARIABLE n.L                   =        4.000  number of open facilties

    ----     57 VARIABLE isopen.L  use facility

    facility1 1.000,    facility2 1.000,    facility3 1.000,    facility4 1.000


    ----     60 PARAMETER locations  

                        x           y

    facility1      26.707      31.796
    facility2      68.739      68.980
    facility3      28.044      67.880
    facility4      76.921      34.929

有关详细信息，请参阅 here。

基本上我们解决了两个模型：

模型 1 找到所需的仓库数量（最小化受最大距离限制的数量）
模型 2 找到仓库的最佳位置（总距离最小）

解决模型 1 后我们看到（对于 50 个客户的随机问题）：我们需要三个仓库。虽然没有 link 超过最大距离限制，但这不是最佳放置。解决模型 2 后，我们看到：这现在通过最小化 link 的长度总和来优化放置三个仓库。准确地说，我最小化了平方长度之和。去掉平方根让我可以使用二次求解器。

两个模型都是凸混合整数二次约束问题类型 (MIQCP)。我使用现成的求解器来求解这些模型。

Answer 2

Python 以 Gurobi 作为求解器的代码：

from gurobipy import *
import numpy as np
import pandas as pd
import networkx as nx
import matplotlib.pyplot as plt



customer_num=15
dc_num=10
minxy=0
maxxy=10
M=maxxy**2
max_dist=3
service_level=0.7
covered_customers=math.ceil(customer_num*service_level)
n=0
customer = np.random.uniform(minxy,maxxy,[customer_num,2])


#Model 1 : Minimize number of warehouses

m = Model()

###Variable
dc={}
x={}
y={}
assign={}

for j in range(dc_num):
    dc[j] = m.addVar(lb=0,ub=1,vtype=GRB.BINARY, name="DC%d" % j)
    x[j]= m.addVar(lb=0, ub=maxxy, vtype=GRB.CONTINUOUS, name="x%d" % j)
    y[j] = m.addVar(lb=0, ub=maxxy, vtype=GRB.CONTINUOUS, name="y%d" % j)

for i in range(len(customer)):
    for j in range(len(dc)):
        assign[(i,j)] = m.addVar(lb=0,ub=1,vtype=GRB.BINARY, name="Cu%d from DC%d" % (i,j))

###Constraint
for i in range(len(customer)):
    for j in range(len(dc)):
        m.addConstr(((customer[i][0] - x[j])*(customer[i][0] - x[j]) +\
                              (customer[i][1] - y[j])*(customer[i][1] - \
                              y[j])) <= max_dist*max_dist + M*(1-assign[(i,j)]))

for i in range(len(customer)):
    m.addConstr(quicksum(assign[(i,j)] for j in range(len(dc))) <= 1)

for i in range(len(customer)):
    for j in range(len(dc)):
        m.addConstr(assign[(i, j)] <= dc[j])

for j in range(dc_num-1):
    m.addConstr(dc[j] >= dc[j+1])

m.addConstr(quicksum(assign[(i,j)] for i in range(len(customer)) for j in range(len(dc))) >= covered_customers)

#sum n
for j in dc:
    n=n+dc[j]

m.setObjective(n,GRB.MINIMIZE)

m.optimize()

print('\nOptimal Solution is: %g' % m.objVal)
for v in m.getVars():
    print('%s %g' % (v.varName, v.x))
#     # print(v)


# #Model 2: Optimal location of warehouses

optimal_n=int(m.objVal)
m2 = Model()   #create Model 2

# m_new = Model()

###Variable
dc={}
x={}
y={}
assign={}
d={}

for j in range(optimal_n):
    x[j]= m2.addVar(lb=0, ub=maxxy, vtype=GRB.CONTINUOUS, name="x%d" % j)
    y[j] = m2.addVar(lb=0, ub=maxxy, vtype=GRB.CONTINUOUS, name="y%d" % j)

for i in range(len(customer)):
    for j in range(optimal_n):
        assign[(i,j)] = m2.addVar(lb=0,ub=1,vtype=GRB.BINARY, name="Cu%d from DC%d" % (i,j))

for i in range(len(customer)):
    for j in range(optimal_n):
        d[(i,j)] = m2.addVar(lb=0,ub=max_dist*max_dist,vtype=GRB.CONTINUOUS, name="d%d,%d" % (i,j))

###Constraint
for i in range(len(customer)):
    for j in range(optimal_n):
        m2.addConstr(((customer[i][0] - x[j])*(customer[i][0] - x[j]) +\
                              (customer[i][1] - y[j])*(customer[i][1] - \
                              y[j])) - M*(1-assign[(i,j)]) <= d[(i,j)])
        m2.addConstr(d[(i,j)] <= max_dist*max_dist)

for i in range(len(customer)):
    m2.addConstr(quicksum(assign[(i,j)] for j in range(optimal_n)) <= 1)

m2.addConstr(quicksum(assign[(i,j)] for i in range(len(customer)) for j in range(optimal_n)) >= covered_customers)

L=0
L = quicksum(d[(i,j)] for i in range(len(customer)) for j in range(optimal_n))

m2.setObjective(L,GRB.MINIMIZE)

m2.optimize()


#########Print Optimization Result
print('\nOptimal Solution is: %g' % m2.objVal)

dc_x=[]
dc_y=[]
i_list=[]
j_list=[]
g_list=[]
d_list=[]
omit_i_list=[]
for v in m2.getVars():
    print('%s %g' % (v.varName, v.x))
    if v.varName.startswith("x"):
        dc_x.append(v.x)
    if v.varName.startswith("y"):
        dc_y.append(v.x)
    if v.varName.startswith("Cu") and v.x == 1:
        print([int(s) for s in re.findall("\d+", v.varName)])
        temp=[int(s) for s in re.findall("\d+", v.varName)]
        i_list.append(temp[0])
        j_list.append(temp[1])
        g_list.append(temp[1]+len(customer))  #new id mapping to j_list
    if v.varName.startswith("Cu") and v.x == 0:
        temp=[int(s) for s in re.findall("\d+", v.varName)]
        omit_i_list.append(temp[0])
    if v.varName.startswith("d") and v.x > 0.00001:
        d_list.append(v.x)


#########Draw Netword
# prepare data
dc_cor=list(zip(dc_x,dc_y))
dc_list=[]
for i,k in enumerate(dc_cor):
    temp=len(customer)+i
    dc_list.append(temp)

df=pd.DataFrame({'Customer':i_list,'DC':j_list,'DC_drawID':g_list,'Sqr_distance':d_list})
df['Sqrt_distance']=np.sqrt(df['Sqr_distance'])
print(df)

dc_customer=[]
for i in dc_list:
    dc_customer.append(df[df['DC_drawID'] == i]['Customer'].tolist())
print('\n', dc_customer)

#draw
G = nx.DiGraph()

d_node=[]
e = []
node = []
o_node = []
for c, k in enumerate(dc_list):
    G.add_node(k, pos=(dc_cor[c][0], dc_cor[c][1]))
    d_node.append(c)
    v = dc_customer[c]
    for n, i in enumerate(v):
        G.add_node(i, pos=(customer[i][0], customer[i][1]))
        u = (k, v[n])
        e.append(u)
        node.append(i)
        G.add_edge(k, v[n])
for m,x in enumerate(omit_i_list):
    G.add_node(x, pos=(customer[x][0], customer[x][1]))
    o_node.append(x)
nx.draw_networkx_nodes(G, dc_cor, nodelist=d_node, with_labels=True, width=2, style='dashed', font_color='w', font_size=10, font_family='sans-serif', node_shape='^',
node_size=400)
nx.draw_networkx_nodes(G, customer, nodelist=o_node, with_labels=True, width=2, style='dashed', font_color='w', font_size=10, font_family='sans-serif', node_color='purple',
node_size=400)
nx.draw(G, nx.get_node_attributes(G, 'pos'), nodelist=node, edgelist=e, with_labels=True,
        width=2, style='dashed', font_color='w', font_size=10, font_family='sans-serif', node_color='purple')


# Create a Pandas Excel writer using XlsxWriter as the engine.
writer = pd.ExcelWriter('Optimization_Result.xlsx', engine='xlsxwriter')

# Convert the dataframe to an XlsxWriter Excel object.
df.to_excel(writer, sheet_name='Sheet1')
writer.save()

plt.axis('on')
plt.show()

设施位置 - 最小化为具有距离限制的客户提供服务的设施的算法

Facility Location - Algorithm to Minimize facilities serving customers with distance constraint

algorithm

optimization

logistics