当前位置：首页 > news >正文

论文复现丨物流中心选址问题：蜘蛛猴算法求解

news 来源：原创 2024/9/21 5:40:14

路径优化系列文章：

1、路径优化历史文章
2、物流中心选址问题论文复现丨改进蜘蛛猴算法求解

物流中心选址问题

一般物流中心选址问题是指：在有限的用户(即需求点)中找出一定数量的地点建立配送中心，实现从物流中心到用户之间的配送，使得各配送中心到由它配送用户的需求量与距离的乘积之和最小。

本文假设：

（１）不考虑物流中心的规模及经济问题；

（２）配送中心必须满足所有位置的需求；

（３）不考虑其他费用。

模型

推文的测试模型包括基准函数和中心选择问题，以下是5个基准函数，可以任意增加：

在文件HS：

基准函数纬度n目前设置为20，可修改。

物流中心选址问题的目标函数和约束见参考文献：

数据

论文数据，31个城市选6个中心的数据（复制pdf数据到txt）：

1 (1304,2312) 20 17 (3918,2179) 90 
2 (3639,1315) 90 18 (4061,2370) 70 
3 (4177,2244) 90 19 (3780,2212) 100 
4 (3712,1399) 60 20 (3676,2578) 50 
5 (3757,1187) 70 21 (4029,2838) 50 
6 (3326,1556) 70 22 (4263,2931) 50
7 (3238,1229) 40 23 (3429,1908) 80 
8 (4196,1044) 90 24 (3507,2376) 70 
9 (4312,790) 90 25 (3394,2643) 80 
10 (4386,570) 70 26 (3439,3201) 40 
11 (3007,1970) 60 27 (2935,3240) 40 
12 (2562,1756) 40 28 (3140,3550) 60 
13 (2788,1491) 40 29 (2545,2357) 70 
14 (2381,1676) 40 30 (2778,2826) 50 
15 (1332,695) 20 31 (2370,2975) 30 
16 (3715,1678) 80

按行和空格拆分数据，具体读取在文件read_decode_paint。

编码解码

编码

基准函数和选址问题都需要生成指定纬度下的随机初始解，2-16生成方式：

基准的最大最小值见函数定义域，对于选址问题，设置最大值为1，最小值为0。

选址问题的初始解还需转换编码：

这其中会出现问题，候选中心数大于或小于设定值。文献是选择6个点，当编码1的计数大于或者小于6，是非可行解，需要修正解。

修正方法：大于6个点时，排序编码x的所有基因值，更正最大6个以外的1为0；小于6个点时，排序编码x的所有基因值，更正不足6个外的0为1.。

如3个点选择2个中心

x=[0.6,0.55,0.8]，转换编码 y=[1,1,1]，修正编码为y= [1,0,1]，保留最大2个，0.55对应的1修正为0

x=[0.35,0.4,0.8]，转换编码 y=[0,0,1]，修正编码为y= [0,1,1]，补充次大的1个，0.4对应的0修正为1

核心代码：

def f_xuan(n,num,x):if x==[]:                             # 为空时生成初始解x = []for i in range(n):xx = np.random.random()       # 0到1生成初始解x.append(xx)y = []for xx in x:                          # 编码转化为0,1if xx < .5:y.append(0)else:y.append(1)idx = np.array(x).argsort()           # 排序xy1 = np.array(y)if sum(y)<num:                        # 如果1的计数小于候选数，补充1y1[idx[-num:]] = 1else:                                 # 否则删除1y1[idx[:-num]] = 0return x, y1.tolist()

解码

基准函数解码较为简单，新解传入对应函数就能返回目标。选址问题设计1种贪婪选择策略：对于每个点，计算点到所有候选中心的距离与需求量的乘积，选择乘积最小对应的中心：

代码如下：

def decode(y,inf):center_idx = [i for i in range(len(y)) if y[i] == 1]  # 候选位置读取center_road = [[] for i in range(len(center_idx))]f = 0for w in range(len(y)):if y[w] == 0:disdemand = [np.sqrt((inf[w][0]-inf[c][0])**2 + (inf[w][1]-inf[c][1])**2)*inf[w][2] for c in center_idx]f += min(disdemand)                          # 目标更新idx = disdemand.index(min(disdemand))        # 最小目标对应的候选中心center_road[idx].append(w)                   # 对应候选中心加入点

改进算法

算法步骤

Step 1. 设置参数。种群规模 P ，最大迭代次数 N ，干扰率 pr ，最大组数 MG ，全局领导限制 gll，局部领导限制lll ，指引因子最大值cmax 及最小值cmin ，变异概率β。

Step 2. 用公式(2-16)生成群体中所有蜘蛛猴的初始位置。然后对每个蜘蛛猴的适应度值进行计算并排序，选出局部和全局领导者。

Step 3. 当不满足终止条件时，执行以下步骤：

Step 3.1. 为了寻找最优解，利用(3-2)式产生新的位置；结合每只猴子适应性的大小择优选出适宜的位置。

Step 3.2. 根据(2-19)式，得出全部组中个体被筛选的可能性；然后利用(3-4)式将选中的蜘蛛猴进行位置更新。

Step 3.3. 对全部组的猴子筛选，选出局部及全局指引者的位置。

Step 3.4. 如果任何一个局部指引者在指定次数(即局部指引限制lll )没有更新其位置，那么局部指引者会再次带领该组全部猴子寻找食物，根据扰动率 pr 的大小采取式(2-21)来更新它们的位置。

Step 3.5. 执行 Cauchy 变异操作，根据式(3-5)更新蜘蛛猴的位置，选出较优的位置。

Step 3.6. 如果全局指引者在指定次数(即全局领导限制 gll )没有更新其位置，且此时没有达到最大组数 MG ，最优指引者会将群体进行分组；否则，最优指引者会把全部小组并为一组。

3-2更新方式：

代码没有加λ，可以自行调整。

3-4更新方式：

3-5更新方式：

核心代码：

for g in range(self.N):if g == 0:Total_x, answer = self.creat(wd, hs)                      # 2-16生成初始位置            best_idx = answer.index(min(answer))gbest, pbest = Total_x[best_idx], Total_x[best_idx]       # 初始全局和局部领导者pbidx = [best_idx]pbcount = [0]result[0].append(g)result[1].append(min(answer))for w in range(self.P):r = np.random.randint(0, self.P)                          # 随机选择索引xr = Total_x[r]x = Total_x[w]mu = 1 - g/self.Nxnew = x + mu*(pbest - x) + mu*(xr - x)                  # 3-2更新方式prob = .9*(answer[w]/max(answer)) + .1              if np.random.random() < prob:                            # 2-19方法概率计算，x = xnewf = answer[w]fmin = min(answer)fave = sum(answer)/len(answer)if f<= fave:c = self.cmin + (self.cmax-self.cmin)/(fave - fmin)*(f - fmin)  # c值计算方式else:c = self.cmaxxnew = x + c*np.random.random()*(gbest - x) + c*np.random.uniform(-1,1)*(xr - x)# 3-4计算方式if hs == 'f_xuan':x, y = rdp.f_xuan(31,6,xnew)center_road,center_idx,f = rdp.decode(y,inf)else:x1, f = eval(f'HS.{hs}({wd},{xnew.tolist()})')if f < answer[w]:                                   # 新解优于原解时更新Total_x[w] = xnewanswer[w] = fbest_idx1 = answer.index(min(answer))if best_idx1 == best_idx:                               # 全局最优位置不变时，gl加1gl += 1else:                                                   # 否则归0gl = 0lenth = int(self.P/len(pbidx))                          # 目前的组数for w in range(len(pbidx)):pbans = answer[lenth*w:lenth*(w+1)]                 # 截取每个组的目标pbid = pbans.index(min(pbans))                      # 每个组的最优解（局部领导）if pbid == pbidx[w]:                                # 局部领导位置不变时，该位置计数加1pbcount[w] += 1else:                                               # 否则归零pbcount[w] = 0if pbcount[w] == self.lll:                          # 如果达到局部领导限制for k in range(lenth):if np.random.random() < pr:                 # 一定概率下随机生成初始解xnew, f = eval(f'HS.{hs}({wd},[])')else:r = np.random.randint(0, P)                         # 否则随机选择索引xr = Total_x[r]          x = Total_x[w]xnew = x + np.random.random()*(pbest - x) + np.random.random()*(xr - x)  # 2-21更新方式f = HS.f1c(xnew)Total_x[lenth*w + k] = xnewanswer[lenth*w + k] = ffor w in range(self.P):                                 #  Cauchy 变异操作x = Total_x[w]if np.random.random() < self.β:a = 2*(1- g/self.N)xnew = x + a*math.tan(np.pi*(np.random.random() - .5))  # 3-5更新方式else:xnew = xif hs == 'f_xuan':x, y = rdp.f_xuan(31,6,xnew)center_road,center_idx,f = rdp.decode(y,inf)else:x1, f = eval(f'HS.{hs}({wd},{xnew.tolist()})')if f < answer[w]:                                   # 新解优于原解时更新Total_x[w] = xnewanswer[w] = fif gl == self.gll:                                   # 如果到达全局领导限制，组数加1            mg += 1if mg == self.MG:                                    # 如果组数到达限制，归1，全局领导和局部领导相同mg = 1best_idx = answer.index(min(answer))pbidx = [best_idx]else:                                               # 没到达组数限制时pbidx = []lenth = int(self.P/mg)                          # 新的组数（+1）对应的组长度pbcount = []for mgg in range(mg):               pbans = answer[lenth*mgg:lenth*(mgg+1)]     # 组长度对应的目标bid = pbans.index(min(pbans))               # 局部领导位置pbidx.append(bid)pbcount.append(0)                           # 局部领导位置计数

结果

代码运行环境
windows系统，python3.6.0,第三方库及版本号如下：

matplotlib (3.3.4)
numpy (1.19.5)
pandas (1.1.5)
xlrd (1.2.0)
XlsxWriter (3.0.3)

第三方库需要在安装完python之后，额外安装，以前文章有讲述过安装第三方库的解决办法。

主函数

设计主函数如下：

from cmasmo import CMASMO
import paint
import HS
import read_decode_paint as rdp
import numpy as npP = 50       # 种群规模       
N = 1000     # 迭代次数
pr = .5      # 扰动率
MG = 5       # 最大组数
gll = 50     # 全局领导限制
lll = 1500   # 局部领导限制
cmax = .8    # 最高学习因子
cmin = .3    # 最低学习因子
β = .8       # 变异概率wd = 20      # 函数纬度，目前是20
hs = 'f1'    # 函数:f1到f5 和 f_xuanif hs == 'f_xuan':  # 选址问题视数据而定，目前是31wd = 31cm = CMASMO(P, N, pr, MG, gll, lll, cmax, cmin, β)
x,result = cm.total(wd, hs)print('*'*20)
if hs == 'f_xuan':inf = rdp.get()x, y = rdp.f_xuan(31,6,x)              # 目前是31点选择6个中心，修正时cmasmo对应行也修正center_road,center_idx,f = rdp.decode(y,inf)print('选择中心编码：',y)print('对应点：',np.array(center_idx)+1)print('目标：',f)rdp.draw(center_road,center_idx,f,inf)
else:x1, f = eval(f'HS.{hs}({wd},{x})')print('函数解：',x)print('目标：',f)
rdp.draw_change(result)

运行结果

f1结果如下：

f1目标随迭代次数的变化图如下：

f2目标随迭代次数的变化图如下：

f3目标随迭代次数的变化图如下：

f4目标随迭代次数的变化图如下：

f5目标随迭代次数的变化图如下：

f_xuan结果如下：

选址图

f_xuan目标随迭代次数的变化图如下：

结论

本推文主要是对参考论文进行复现，为了解可行设计了选址中心的编码修正方式，为了求解效率设计了贪婪的中心选择方式。其中基准函数容易陷入局部最优，可从算法参数、算法流程、初始解或者函数定义域等进行优化，基准函数目前是5个，其余的可自行在HS添加。
数据和算法皆可修改（修改选址数据时，其一行包含1个或者2个点的相关数据，否则报错，中心点和选择个数在main和cmasmo都需修改）。
参考文献：蜘蛛猴算法的改进及其在物流中心选址问题中的应用_杨转（第3章）