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

Question

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

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

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

score 5 · Accepted Answer

这属于设施选址问题的范畴。关于这些问题的文献非常丰富。p 中心问题与您想要的很接近。

一些注意事项：

除了求解正式的数学优化模型外，还经常使用启发式（和元启发式）。
距离是实际旅行时间的粗略近似值。这也意味着近似解决方案可能已经足够好了。
除了找到服务所有客户所需的最少设施数量外，我们还可以通过最小化距离来优化位置。
纯“最小化设施数量”的数学规划模型可以表述为混合整数二次约束问题 (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

有关更多详细信息，请参见此处。

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

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

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

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

score 3 · Accepted Answer

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

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

2 回答 2

Related

Reference