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

Question

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

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

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

Answer 1

这属于设施位置问题。关于这些问题的文献非常丰富。 p-center 问题与您想要的很接近。

一些注意事项：

1. 除了求解正式的数学优化模型外，还经常使用启发式（和元启发式）。
2. 距离是实际行程时间的粗略近似值。这也意味着近似解可能就足够了。
3. 除了找到服务所有客户所需的最少设施数量外，我们还可以通过最小化距离来优化位置。
4. 纯“设施数量最小化”的数学编程模型可以表述为混合整数二次约束问题 (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. 模型 1 找到所需的仓库数量（根据最大距离限制最小化数量）
2. 模型 2 找到仓库的最佳位置（最小化总距离）

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

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

Answer 2

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

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()