合理的仓库选址对现代物流管理有着重要影响，通过对仓库进行合理规划，不仅可以节约配送成本，还可以提高客户满意度及企业的核心竞争力[1]。仓库选址问题是一个经典的组合优化问题，该问题通过确定设施点的位置，将需求点分配给开放的设施点，确保每个需求点的需求均能被满足，并使总成本最小化。现如今，有容量的设施选址问题已成为物流管理、供应链管理及运筹优化领域的研究热点。李小川等[2]设计了改进的烟花算法对物流配送中心选址问题进行求解；李捷承等[3]针对物流配送设施选址的特点，设计了一个基于BIRCH(Balanced Iterative Reducing and Clustering Using Hierarchies）聚类的物流设施选址算法，可较大程度地节约长期运营成本；黄凯明等[4]对于物流网络多层级设施选址—路径规划问题，建立了混合整数规划数学模型，并提出量子进化算法与遗传算法协同的双智能算法集成求解方案；陈勇等[5]建立了逆向物流网络模型，考虑废旧家电回收的数量、质量、客户需求量、回收中心对居民产生的负效用，为企业在物流设施选址及不同周期下的市场缺货量、设施间流量分配等提供灵活的决策方案；张雨等[6]以时间路径优化为研究对象，对聚类算法和最小支撑树法进行改进，构建物流中心选址和物流配送区域划分模型，以提升物流企业服务效率与服务质量。在物流仓库选址模型中，考虑年运行成本及仓库容量约束的研究较少。鉴于此，文中建立的有容量的仓库选址模型考虑了仓库年固定成本、车辆行驶成本和仓库容量约束，使得模型更具有实际应用意义。

P-中位选址问题（P-Median）是一类典型的组合优化问题，属于NP-hard问题，在设施选址、交通、物流等领域有着较为广泛的应用[7,8]。本文所建立的有容量的应急医疗设施选址模型是对P-中位选址问题的一种扩展，也属于NP-hard问题。对于NP-hard问题的求解，精确算法很难在可接受的计算时间范围内提供最优解[9]。虽然大多数元启发式算法不能保证得到最优解，但它们可以在合理的时间内提供良好结果。已经有不少学者将遗传算法[10,11,12]、禁忌搜索算法[13]、蚁群优化算法[14]、模拟退火算法[15]、免疫优化算法[16]等元启发式算法应用到选址问题中。免疫优化算法被发现是一种解决多种优化问题的高效算法[17]，其利用免疫机制保持种群多样性。此外，免疫优化算法还具有自组织和鲁棒性强的优点，适用于求解选址优化问题[18,19,20]。但是，免疫优化算法的免疫操作较为简单，在求解过程中容易陷入局部最优。因此，根据仓库选址问题的特点，对免疫优化算法设计两种交叉算子与变异算子，使得算法在迭代过程中平衡局部搜索和全局搜索。此外，对算法中的参数进行测试，选取最优参数组合方式，提高算法性能，降低算法求解运行时间。

本文建立的仓库选址模型考虑了设施点的容量、年建设成本和车辆行驶成本。有容量的仓库选址模型需从候选设施点集中选取部分作为设施点，为需求点提供货物运输服务。保证每个需求点均由设施点提供服务，且设施点服务的总需求量不超过设施点的服务容量限制。在选址问题中，车辆的行驶费用与设施点的建造费用往往是相互矛盾的。建立较少的仓库可节省较多的建设费用，但会增加车辆的交通成本，反之亦然。因此，构建适当数量的仓库可以帮助平衡仓库建设成本与车辆行驶成本。

为便于选址模型建立，给出以下假设：(1)一个设施点可为多个需求点提供服务，但一个需求点仅由一个设施点提供服务；(2)车辆行驶成本与行驶距离之间存在简单线性关系；(3)区域内的需求点拟用离散型变量分布点集N表示，候选设施点从点集M中产生；(4)设施点仓库建设为统一标准。

模型所用符号说明如下：

集合：

N：需求点的集合；

M：设施点的集合；

参数：

dij:i需求点到j设施点之间的距离，i∈N,j∈M;

p：设施点的个数；

B：设施点的建设成本；

r：计算设施点的折现率；

t：新建设一个设施点的使用年限；

k：道路曲折系数；

Di：需求点i的需求量；

G：设施点最大可服务的容量。

决策变量：

基于上述假设、定义的参数和决策变量，将仓库选址问题的数学模型表示为：

约束条件：

目标函数（1）为设施点的年建设成本与车辆行驶成本之和。固定资产的价值可能会随着时间而变化。在仓库的使用年限内，将固定成本折算为年成本更为合理。贴现率r常用来表示建筑成本的当前价值与未来价值的关系[21,22]。因此，目标函数的第一部分为设施点的年建设成本，由建设成本与系数相乘得到[21]。目标函数的第二部分是车辆行驶成本。约束条件（2）表示仅当设施点j开放时，设施点j才可为需求点i提供运货服务。约束条件（3）表示需求点i仅由一个设施点提供服务。约束条件（4）表示设施点开放的个数为p。约束条件（5）为服务容量的约束，即设施点所服务的容量不可超过其容量上限。约束条件（6）、（7）是对决策变量xij、yj的取值范围进行约束。

免疫优化算法是一种基于生物免疫系统的智能优化算法。在该算法中，模型中的目标函数和解分别对应于生物免疫系统中的“抗原”和“抗体”[18]。在算法迭代过程中，抗体与抗原之间的亲和值逐渐提高，最终生成亲和值最优的抗体。免疫优化算法流程如图1所示。免疫优化算法在免疫选择的作用下，反映了免疫记忆功能和免疫系统的自我调节功能。

图1 免疫优化算法流程   

Fig.1 Flow chart of immune optimization algorithm

考虑到有容量仓库选址模型的具体特征，设计改进的免疫优化算法，使其更适合有容量的仓库选址模型，并克服经典免疫算法易陷入局部最优解的问题。改进的免疫优化算法主要由3部分构成：(1)初始化：确定解的表达形式，随机产生初始抗体，将满足条件的一定数量的抗体保存，以便算法能够快速地寻优迭代；(2)抗体评价：抗体评价是通过计算抗体的亲和力评价抗体优劣，亲和力计算包括抗体与抗原的亲和值、抗体密度。将质量高的若干抗体存入记忆库中，以便算法快速收敛到最优解；(3)遗传操作：在遗传操作中，利用单点交叉、二点交叉、变异算子生成子代抗体，既加强了对有希望空间的搜索利用，又能保持种群的多样性，使其快速收敛到最优解。

仓库选址模型的选址方案可以形成长度为l的抗体（l为开放设施点的数量），每个抗体表示开放设施点的序列，即解的序列。对于m个候选设施点中选取l个开放设施的选址问题，假设m个候选设施点的编号分别为1,2,…,m (m>1)，则抗体l=[M1,M2,…,Ml]表示一个可行的解决方案，其中，M1,M2,…,Ml表示1,2,…,m中不重复的l个序号。

与经典免疫优化算法相同，初始解首先通过洗牌方式产生：将候选设施点集1,2,…,m (m>1)的顺序打乱，随机选择两个不同位置的数字并进行交换，执行次。将交换后序列的前l个数字作为初始抗体。但由此产生的抗体具有随机性，抗体质量较差。因此，本文对抗体进行判断，选择优质抗体作为初始抗体。抗体中的每个设施点为Mj,(j∈M),usdcap(j)为设施点j已服务的需求量，demandi为需求点i的需求量，Cj为j设施点的服务容量。将需求点i分配给离其最近且满足容量约束（usdcap(j)+demandi≤Cj）的设施点j，令xij=1，并更新usdcap(j)=usdcap(j)+demandi；遍历完所有需求点i，如果demandi均能被满足，且设施点所服务的总需求不超过该设施点的容量限制，则判定抗体合格；否则，重新产生新的抗体并检验是否合格。本文在算法中设置记忆库存储高质量的抗体，以保留优秀抗体的优秀特征，初始记忆库中的抗体是随机生成的合格抗体。将合格的P1个抗体与记忆库中P2个抗体组成初始抗体群。

抗体与抗原之间的亲和值表示抗体对抗原的识别程度。通过式（8）计算抗体与抗原间的亲和值，即适应度值。

抗体间的相似值表示它们之间的相似性。本文采用Smith等[23]提出的R-bit连续法计算抗体相似值。当两个抗体有连续R位或更多位置上具有相同编码时，表示两种抗体相似，否则表示两种抗体不相似。

其中，tk,v为抗体k、v之间的相同位数，l为每个抗体的长度。例如，抗体xk=[53 23 12 45 61 20 37 56 24 65 34 1]和xv=[25 39 53 32 11 47 12 16 15 19 45 5]的相同位数为3，则xk和xv之间的Sk,v为0.25。

抗体密度φk，即抗体浓度，表示群体中相似抗体的比例。

如果Sk,v>R，则Bk,v=1，否则Bk,v=0。R是预定义的阈值[23]。

通过式（11）可计算抗体的期望繁殖概率，式（11）由抗体的亲和值和抗体密度两部分构成。

其中，η为[0,1]区间内的常数。

由式（11）可知，抗体与抗原的亲合值越大，期望繁殖概率越高；抗体密度越大，期望繁殖概率越小。因此，式（11）不仅可以鼓励选择适应度值好的抗体，而且可以抑制抗体间的密度，从而确保种群的多样性。当免疫优化算法抑制高密度的抗体时，部分高密度的抗体可能已接近于最优解，这样会导致已求得的最优解丢失。因此，本文采用精英保留策略：在更新记忆库时，先将与抗原亲和值高的P个抗体存入记忆库中，保留抗体群中精英个体的优秀特征；再提取前P1个抗体作为父代群体，按照期望繁殖概率进行选择和遗传操作，以保持优良抗体的特性。

通过轮盘赌的机制进行选择操作，轮盘选择法是一种比较受欢迎的选择方法。根据抗体的期望繁殖概率，在父代群体P1中选择两个抗体。期望繁殖概率越大，被选择的概率就越高。

交叉在遗传操作中起着至关重要的作用，它可以提高抗体间的差异度，以便算法快速收敛到最优解。考虑是在整数编码的遗传操作环境下，本文使用时间复杂度低的两种交叉操作方式：单点交叉和两点交叉。对于两种交叉方式的参数选择，设置交叉率的范围为θ∈[θmin,θmax]，随机生成[0,1]之间的随机数γ，若γ≤θ，则选择单点交叉，反之，选择两点交叉操作。

单点交叉：给定两个抗体parent1、parent2，首先在2和l-1之间随机选择一个交叉点q。交换两个抗体在q点之后的元素，得到parent1,、parent2,。若parent1,中更改的序号r在parent1中已存在，则找到parent1中序号r所在的位置R，并找到parent2中R位置的序号r,，将parent1,中的序号r改为r,。重复此步骤，直到抗体中无重复序号，得到子代抗体child1和child2。

单点交叉操作如图2所示。对于亲本抗体parent1=[102 21 13 9 17 6 23 4]和parent2=[11 8 20 1 15 12 21 14 5]，随机选择交叉点q=6。将两个抗体第6个位置后的元素交换得到parent1,=[10 2 21 13 9 12 21 14 5],parent2,=[11 8 20 115 17 6 23 4]。在parent1,中有重复元素21，因此找到parent1中元素21所在的第3个位置，并找到parent2中位置3的序号为20，将parent1,第7个位置的元素改为20。最后，得到无重复元素的子代抗体child1=[10 2 21 13 9 12 20 145]和child2=[11 8 20 1 15 17 6 23 4]。

图2 单点交叉例子   

Fig.2 An example of one-point crossover

两点交叉：给定两个亲本抗体parent1、parent2，首先在2和l-1之间随机选择两个不同的交叉点q1、q2，交换两个抗体交换q1和q2之间的元素，得到抗体parent1,、parent2,。与单点交叉相同，需要消除子代抗体中的重复元素，得到子代抗体child1、child2。

两点交叉操作如图3所示。对于亲本抗体parent1=[102 21 13 9 17 6 23 4]和parent2=[11 8 20 1 15 12 21 14 5]，随机选择两个交叉点q1=4,q2=6。将两个抗体q1和q2之间的元素交换，得到parent1,=[10 2 21 1 5 12 6 23 4],parent2,=[11 8 20 13 8 17 19 14 5]。parent2,中有重复元素8，因此找到parent2中元素8所在的第2个位置，并找到parent1中第2个位置的元素为2，将parent2,第5个位置的元素改为2。最后，得到无重复元素的子代child1=[10 2 21 1 5 12 6 23 4]和child2=[11 8 20 13 2 17 19 14 5]。

突变指通过在抗体中引入随机变异以增加种群多样性，消除算法在无希望区域的停滞，探索新搜索区域的过程[24]。对于一个抗体，将其含有的元素记为child，即开放设施点的集合；抗体中不包含的元素记为Unopen，即未开放的设施点集合。随机选择候选设施点u∈child，未开放的设施点uˉ∈Unopen，将抗体中的u改为uˉ以实现突变操作。为了增加种群多样性，重复此过程，直到突变率（突变次数）达到μ。

在实现突变的过程中，还应考虑保护优秀解的特征，因为在迭代过程中，可能已求得一些解接近于最优解，直接对解执行突变操作，可能会丢失已求得的优秀解[24]。因此，在突变过程中选择适应度值fitness最好的进行下一次迭代。如图4所示，假设抗体长度l=3，突变率μ=3，当候选设施点集为M={1,2,…,12}时，则需要对抗体执行3次突变操作，记录每次突变操作的适应度值，选择适应度值最优的抗体进入下一次迭代。

图3 两点交叉的例子   下载原图

Fig.3 An example of two-point crossover

Fig.4 An illustration of the mutation withμ=3,p=3 and N={1,2,⋯,12}   

图4 μ=3,l=3,N={1,2,…,12}的突变操作例子

以上海市S物流公司为案例，根据实地调研，S物流公司每天由仓库向各网点运输货物。该公司在上海市共有51个网点分布，因此获取51个需求点的坐标与需求量，并选取15个候选仓库点，需求点的需求量为每天运输的车辆数。模型中各参数如表3所示。需求点与设施点的相关数据如表1、表2所示。

根据经纬坐标计算两点之间的距离，使用式（12）将经纬距离转换为两个节点i和j之间的实际行驶距离dij:

其中，(xi,yi)、(xj,yj)为两个节点的经度和纬度坐标，6 370为地球半径（km）。通过式计算出两点间的线性距离。抽取30组两点间的线性距离，将其与百度地图得到的实际行驶距离相比较，得到误差值k。

参数的合理设置对算法的有效性和计算效率有着重要影响，算法中的参数：种群规模P=P1+P2(P1为父代群体数量,P2记忆库容量）、迭代次数iter、记忆库容量P2、交叉率的区间[θmin,θmax]、期望繁殖概率评估参数η、突变率μ需要进行合理设置。参数调优和参数控制是元启发式算法中确定参数值的两种常用方法[25]。参数调优是将其他参数固定，仅变化需要测试的参数，找到求解效果最好的参数值；参数控制是将几种参数的值进行变化组合，找到最优组合方式[25]。本文使用这两种方法确定各参数值。在参数测试环节中，考虑到仅对一个案例进行测试具有偶然性，因此将案例生成4个不同规模的子案例并分别求解，测试参数对算法性能的影响，以寻找最优参数组合。

迭代次数在算法中起重要作用，迭代次数设置较高可能会浪费算法运行时间，迭代次数设置较低可能会提前结束搜索过程，无法找到最优解。首先使用精确软件CPLEX分别对4个案例进行求解，并得到最优值；然后使用本文设计的免疫优化算法对4个案例分别独立运行10次，取得平均结果与平均运行时间。表4给出了4个案例的规模与求解情况，包括每个案例变量个数、约束条件个数、最优解（Best）、结果平均值（BOV）、10次求解中找到最优解的数量（Found）、结果与最优值的偏离程度（APD=((BOV-best) best)×100%）、运行平均时间（Time(s)）、求得最优解的迭代次数（Iters）。

表1 需求点坐标与需求量 
Table 1 Demands and coordinates of each demand node    

表2 设施点坐标 
Table 2 Coordinates of each facility node     

表3 模型中的其他参数 
Table 3 Other parameters in the model    

由表4结果可知，4个案例在150次迭代内均求得了最优或接近最优解。因此，本文设置迭代次数Iters=150。

设置种群规模，取6个不同的参数值进行测试。表5给出了不同的种群规模下，4个案例进行10次运行的平均测试结果，包括种群规模、平均运行时间和APD（%）。从表5可以看出，当种群规模大于等于30时，运行时间和解的质量较好。随着种群规模的增大，算法计算时间也越长。考虑到解的质量和计算时间，本文设置种群规模P=30。

表4 4个案例运行10次得到的测试结果 
Table 4 Test results from 10 runs in 4 cases     

表5 种群规模测试结果 
Table 5 Test results of population size     

记忆库用于存储每次迭代产生的精英抗体。精英抗体的数量在种群中占有一定比例。记忆库容量设置过大时，算法容易陷入局部最优；记忆库容量设置过小时，寻找最优解的计算时间会增加。因此，设置合理的记忆库容量，对保留精英抗体的特征与维持种群的多样性具有重要意义。通过设置6个不同的比例，确定最合适的记忆库容量。表6给出了不同的记忆库容量值，4个案例运行10次的平均测试结果，其中，[X]表示不大于X的最大整数。由表6可以看出，当记忆库容量为[0.35*(P）]时，可以得到较好的解决方案。因此，设置P2=[0.35*(P)]，即P=30时，P2=10。

表6 记忆库容量测试结果 
Table 6 Test results of memory capacity    

交叉率用于决定抗体的交叉操作为两点交叉或一点交叉，这对决定如何保留亲本特征较为重要。在算法运行时，与固定交叉率相比，可变的交叉率可以有更大的灵活性去探索更有前景的搜索区域。通过设置18个不同的取值区间，确定最合适的交叉率区间。表7给出了不同的交叉区间下，4个案例运行10次的平均测试结果。由表7所示的测试结果可以看出，在区间[0.5,1]和[0,0.9]内，计算得到的APD（%）值和计算时间较好。

表7 交叉率测试结果 
Table 7 Test results of crossover rate    

变异是遗传操作的主要算子之一，通过变异可以搜索新的区域、逃离当前局部最优情况。确定合适的突变率可以决定抗体变异强度。突变率和交叉率一样重要，会对改进免疫优化算法的计算效率产生较大影响。设置较小的突变率容易在寻优过程中陷入局部最优，设置较大的突变率会失去已得到的优秀特征。表8给出了4个案例在不同突变率下运行10次的平均测试结果。从表8可以看出，当突变率为0.5和0.6时，免疫优化算法的性能更好，可以在较短时间内求得较好的解。

表8 突变率测试结果 
Table 8 Test results of mutation rate     

期望繁殖概率表示对抗体进行评估，参数η用于为抗体对抗原的识别程度、抗体密度赋予不同权重，以此得到抗体的期望繁殖概率。如表9所示，当参数η设置为0.8或0.85时，算法运行时间较短，且APD（%）较低，求得解的质量较高。

表9 参数η测试结果 
Table 9 Test results of parameterη    

不同的参数组合对算法性能也有很大影响。如上所述，交叉率、突变率和参数η有多个较优值，将不同的参数进行组合，得到8种组合方式。使用每种组合方式对4个案例各进行10次运行求解，得到平均运行时间与APD（%）值。不同参数组合的测试结果如表10所示，从测试的运行时间与APD（%）值可以看出，当参数组合[ϑmin,ϑmax]=[0,0.9]、μ=0.5、η=0.8时，算法求解时间与解的质量最优。

表1 0 不同参数组合的测试结果 
Table 10 Test results of different parameter combinations     

基于上述讨论和测试，算法参数设置情况为：迭代次数Iters=150、种群规模P=30、记忆库容量P2=10、交叉率区间[ϑmin,ϑmax]=[0,0.9]、突变率μ=0.5、抗体评估参数η=0.8。

通过上述改进的免疫优化算法及参数设置，对S物流公司的案例进行求解分析。实验环境为Intel(R)Core i5-8250u CPU@1.60GHz 8.00GB内存，操作系统为64位Windows 10，使用Matlab R2017b进行编程。求解得到目标函数为9.116 06×105，运行时间为45.70s。所选择的设施点为3、5、6、8、9、10、11、12、13、15，需求点与设施点分布情况如图5所示（彩图扫OSID码可见）。图5中的黄色圆点为需求点，红色星星为所选的设施点，绿色圆点表示未被选择的设施点，连线表明需求点与设施点之间的关系。从图5可以看出，每个需求点均有指定的设施点进行服务，且需求点均被分配给距离较近的设施点。表11为设施点所服务的需求点和服务的总需求量，可以看出，每个设施点服务的总需求量均不超过服务容量限制。因此，本文的仓库选址模型有效可行，能够解决物流公司的仓库选址与需求分配问题，并有效降低车辆成本与固定成本，提高企业对车辆及仓库的控制。

表1 1 设施点与需求点分配情况 
Table 11 Relationship between facility points and demand points     下载原表

图5 设施点与需求点分布   

Fig.5 Distribution of facility and demand points

为了验证改进免疫优化算法的有效性，本文使用改进免疫优化算法的求解结果与软件CPLEX精确求得的精确解进行对比分析。如表12所示，免疫优化算法与CPLEX求解得到的结果相同，且免疫优化算法运行时间更短。因此，免疫优化算法可有效地解决仓库选址问题。

表1 2 改进的免疫优化算法与CPLEX求解对比 
Table 12 Comparison between improved immune optimization algorithm and CPLEX     

为了验证改进免疫优化算法的求解性能，将其优化结果与经典免疫优化算法的求解结果进行对比分析。图6和图7分别为改进的免疫优化算法与经典免疫优化算法的收敛曲线。由图6和图7可以看出，改进的免疫优化算法无论是收敛速度，还是结果准确度都优于经典免疫优化算法。

图6 经典免疫优化算法收敛曲线   

Fig.6 Convergence curve of classical immune algorithm

图7 改进免疫优化算法收敛曲线   

Fig.7 Convergence curve of improved immune algorithm

本文结合仓库设施选址特点，建立了考虑容量约束的仓库选址模型。根据选址模型特点，设计了改进的免疫优化算法对案例进行求解。在算法中设计了两种交叉算子和变异算子，使得算法可快速收敛到最优解，降低算法运行时间。本文还对改进免疫优化算法的参数进行详细分析，找到最优参数组合方式，提高算法性能。最后，对案例进行求解，验证了模型和算法的有效性与可行性。本文提出的仓库设施选址方法可为决策人员选择最佳布局方案提供科学合理参考。

在下一步研究中，将进一步细化影响仓库选址的各因素，并考虑服务时间窗的约束，更合理地在模型中体现车辆行驶实际情况，以增加模型求解的准确性和实际应用价值。