做物流的人都知道仓库和配送中心的位置选得好不好直接决定运输成本能不能压下来。而选址这件事落到数学上就是一个标准的组合优化问题。我前段时间在MATLAB里用粒子群优化算法跑了一个配送中心选点仿真把几十个需求点、几个候选位置、选几个中心的逻辑完整实现了一遍。这篇文章就把整个思路、代码和踩过的坑都整理出来想用粒子群优化算法做物流选址的朋友可以直接参考。1. 选址问题的本质不是地图上画圈而是组合优化1.1 从一个具体场景说起假设你在负责某个区域的配送网络规划手头有30个客户需求点每个点的位置和需求量是已知的。同时之前整理出6个备选仓库位置但预算有限只能从中选3个建仓。现在要回答一个问题这3个仓选在哪几个备选点上能让“把所有需求点的货物送到对应仓库”的总运输成本最低很多人第一反应是“看地图凭经验”但真实业务里需求点可能有几十上百个候选点错开几条街、偏几百米对整个网络的配送成本影响都非常大。手算几乎不可能所以要用数学模型来描述它。这个问题可以抽象为从候选集合中选出若干个设施点让所有需求点到最近设施点的加权距离总和最小同时把建设成本也算进去。这就是运筹学里常说的设施选址问题。1.2 常见选址模型的取舍用一个数学模型去描述选址主要有几种思路重心法适合单设施选址直接求一个连续坐标使总运输成本最小。优点是简单缺点是现实中你不可能在一个湖中央或者某栋楼的楼顶建仓最后还得人为修正到可用地块上。覆盖模型给每个设施设定一个服务半径看能不能覆盖所有需求点常用于应急中心和消防站。它关心“覆盖没覆盖”不关心“距离多远”对配送成本的计算不够精细。p-中位模型从若干个候选点中选p个设施点使所有需求点到最近设施点的加权距离总和最小。这就是配送中心选址最常用的模型也是本文要处理的问题。p-中心模型选p个设施点让所有需求点到最近设施点的最大距离最小侧重兜底保障适合时效要求严格的场景。p-中位模型的目标函数清晰约束条件也相对简单每个需求点被分配给最近的设施点设施点被选中才承担分配。这个模型看着简单但当候选点数量和需要选择的点数增大时组合爆炸会非常快。6个候选点选3个是20种组合如果变成60个候选点选10个组合数量就变成了天文数字穷举根本跑不动。所以需要启发式算法粒子群就是其中非常实用的一个。1.3 为什么选粒子群和遗传算法、模拟退火的对比我最早接触选址教程时网上清一色用遗传算法。遗传算法确实能解但它的编码、选择、交叉、变异一整套操作下来参数太多调参周期长。模拟退火也不错但对初始温度、降温速率敏感收敛速度偏慢。粒子群优化算法的特点是参数少、编码直观、收敛快特别适合在MATLAB里快速迭代验证。算法核心机制参数数量收敛速度实现难度遗传算法选择、交叉、变异较多中等较复杂模拟退火温度控制、邻域搜索中等偏慢中等粒子群速度更新、位置更新较少较快简单粒子群的另一个优势是连续变量适应性好而选址问题如果把候选点编号看成一个离散整数域也能通过取整映射的方式套到PSO框架里。这就是我这篇文章要展示的核心思路。2. 粒子群优化的运行机制鸟群怎么变成寻优工具2.1 粒子、位置、速度、适应度四个概念的通俗理解粒子群优化算法脱胎于对鸟群觅食行为的模拟。鸟群在天空飞的时候每只鸟既保留自己找到过的最优位置也参考整个鸟群当前发现的最优位置从而调整飞行的方向和速度。放到优化问题里粒子一个粒子就是优化问题的一个候选解。位置候选解中所有决策变量的具体取值。速度下一次迭代中位置的变化方向和幅度。适应度把候选解代入目标函数后得到的评价指标比如总成本。算法每次迭代都更新每个粒子的速度和位置把“自己飞出自己的历史最优位置”和“整体飞向群体最优位置”这两个趋势叠加起来。数学公式如下[ v_{i}^{t1} w v_{i}^{t} c_1 r_1 (pbest_i - x_i^t) c_2 r_2 (gbest - x_i^t) ][ x_{i}^{t1} x_i^t v_{i}^{t1} ]其中 (w) 是惯性权重控制上一代速度对当前速度的影响(c_1)、(c_2) 是学习因子分别控制粒子向自身历史最优 (pbest) 和群体全局最优 (gbest) 学习的强度(r_1)、(r_2) 是0到1之间的随机数。可以这样理解第一项是“惯性”保持原来的飞行趋势第二项是“自我认知”飞向自己曾经找到过的最好位置第三项是“社会认知”飞向整个群体目前发现的最佳位置。三者平衡得当粒子群就会像一群互相交流的鸟逐步逼近最优解。2.2 选址问题如何编码成粒子粒子群的原始形式是为连续优化设计的但选址问题往往要从离散候选点中做选择。常见编码方式有三种方法A连续坐标编码。直接让粒子代表设施点的x,y坐标适合没有候选点限制的重心式选址。缺点是选出来的点可能在实际中无法落地。方法B离散索引映射。本文采用如果有m个候选点要选k个那么粒子就是一个k维连续向量每一维的取值映射到候选点的编号。例如候选点编号从1到6粒子的某个维度取值为4.2向下取整到4就代表选中4号候选点。这种编码实现简单更新公式完全沿用连续PSO只需在计算适应度前做一次取整和去重。方法C二进制编码。每个粒子是一个长度为m的向量每一位0或1表示对应候选点是否被选中。速度通过Sigmoid函数映射到[0,1]区间再转换成概率。这种方法适合候选点数量很大的情况但需要额外处理“选中数量恰好等于k”的约束实现稍显繁琐。我在这个例子里使用方法B原因很简单候选点只有6个选3个用连续索引映射最容易理解而且适应性函数里可以直接用MATLAB的向量运算。2.3 适应度函数的设计运输成本 固定成本的权衡适应度函数是粒子群优化的“指挥棒”。对配送中心选址问题目标函数需要包含两部分运输成本每个需求点分配给离它最近的中心距离乘以需求量然后求和。固定成本某个候选点被选中后需要支付建设或运营成本。一个候选点是否被选中直接决定这部分成本是否计入。数学模型可以写成[ \min ; \sum_{i1}^{N} w_i \times \min_{j \in S} d(i,j) \sum_{j \in S} f_j ]其中 (N) 是需求点数量(w_i) 是需求点 (i) 的需求权重(S) 是被选中的设施点集合(d(i,j)) 是需求点到候选点的距离(f_j) 是候选点 (j) 的固定成本。如果考虑容量限制比如每个中心最多只能服务多少需求那就要在目标函数里加上惩罚项。比如某个中心被分配的总需求超过容量就加上一个很大的惩罚值。这种罚函数法虽然简单但惩罚系数不好标定太小约束不住太大又会让算法搜索困难。后文先不引入容量给一个相对干净的版本。3. 基于MATLAB的粒子群选址实现从建模到出图3.1 算例设定与数据生成为了让代码可以直接跑我设计了一个模拟算例30个随机需求点6个随机候选点从6个候选点中选3个作为配送中心。需求坐标分布在0到100的平面区域内需求量权重在1到10之间候选点固定成本在20到50之间。为了复现方便我用固定随机种子。% 物流配送中心选址 - 粒子群优化算法主脚本 clear; clc; rng(42); % 基本参数 numDemand 30; % 需求点数量 numCand 6; % 候选中心数量 numSelect 3; % 需要选几个中心 % 随机生成需求点坐标x, y 在 [0,100]和需求量权重1~10 demandCoord rand(numDemand, 2) * 100; demandWeight randi([1, 10], numDemand, 1); % 随机生成候选中心坐标 candCoord rand(numCand, 2) * 100; % 候选中心固定建设/运营成本 fixedCost randi([20, 50], numCand, 1); % 粒子群参数 popSize 50; % 种群规模 maxIter 100; % 最大迭代次数 wStart 0.9; % 惯性权重起始值 wEnd 0.4; % 惯性权重结束值 c1 1.5; % 个体学习因子 c2 1.5; % 全局学习因子 % 初始化粒子位置和速度 % 位置范围 (1, numCand1)floor后得到 1~numCand pos rand(popSize, numSelect) * numCand 1; vel rand(popSize, numSelect) * 0.2 - 0.1; % 初始小速度 vmax 0.5 * (numCand - 1); % 速度限幅这里有个细节位置用了rand * numCand 1范围是(1, numCand1)。floor(1.5)1floor(numCand)numCand刚好能覆盖所有候选点编号。如果一开始就用rand(popSize,numSelect)*(numCand-1)1范围是(1, numCand)那floor之后最多只能取到numCand-1最后一个候选点永远不会被选中。这个坑非常隐蔽我第一次跑的时候就是这里出了问题换了随机种子后发现结果总是回避某个候选点排查了半天才找到原因。3.2 适应度计算函数适应度函数要接收“选中的候选点编号集合”然后计算总成本。注意这里的输入是整数索引不是粒子原本的连续值。function cost calcCost(selected, demandCoord, demandWeight, candCoord, fixedCost) % selected: 被选中候选点的编号向量例如 [3 5 6] nDemand size(demandCoord, 1); totalDistCost 0; for i 1:nDemand % 候选点与当前需求点的所有距离 distVector sqrt((demandCoord(i,1) - candCoord(selected,1)).^2 ... (demandCoord(i,2) - candCoord(selected,2)).^2); totalDistCost totalDistCost demandWeight(i) * min(distVector); end totalCost totalDistCost sum(fixedCost(selected)); cost totalCost; end循环里先计算当前需求点到所有已选中心的距离取最小值作为分配距离乘以需求量后累加。这种“最近分配”逻辑非常直观也符合配送中心选址的常规假设每个需求点总是去找离自己最近的中心。固定成本只加一次不管这个中心服务了多少需求点。3.3 离散化处理取整、去重、补齐粒子群更新产生的是连续值但选址需要的是整数编号。我写了一个离散化函数负责把连续向量变成合法的候选点编号集合function sel discreteSelect(pos, numCand, numSelect) % 将连续粒子位置映射为不重复的候选点编号集合 sel floor(pos); % 向下取整 sel max(min(sel, numCand), 1); % 边界修正 selUnique unique(sel); % 去重 if numel(selUnique) numSelect % 去重后数量不足从缺失编号中随机补足 missing setdiff(1:numCand, selUnique); addCount numSelect - numel(selUnique); addIdx missing(randperm(numel(missing), addCount)); sel [selUnique, addIdx]; else % 去重后数量足够只保留前numSelect个 sel selUnique(1:numSelect); end end这个函数的关键在于粒子中的每个维度都可能指向同一个候选点比如某个粒子是 [2.1, 2.9, 5.2]floor之后变成[2, 2, 5]2号候选点重复了。如果不处理后面适应度计算时selected里有重复元素会重复计算固定成本而且可能导致实际选中的点数不足3个。我的做法是先用unique去掉重复再从没被选中的候选点里随机补足。这相当于给候选解加了一个“微扰”虽然随机性会降低局部搜索效率但实现简单而且种群规模足够大时影响不大。3.4 主循环粒子群迭代更新主循环需要保持标准PSO的思想更新速度、更新位置、边界处理、计算适应度、更新个体最优和全局最优。惯性权重采用线性递减让算法前期多探索后期多收敛。% 初始化个体最优和全局最优 pbest pos; pbestCost inf(popSize, 1); for p 1:popSize sel discreteSelect(pos(p,:), numCand, numSelect); pbestCost(p) calcCost(sel, demandCoord, demandWeight, candCoord, fixedCost); end [gbestCost, gbestIdx] min(pbestCost); gbest pos(gbestIdx, :); bestCostHistory zeros(maxIter, 1); for iter 1:maxIter w wStart - (wStart - wEnd) * iter / maxIter; % 线性递减 for p 1:popSize r1 rand(1, numSelect); r2 rand(1, numSelect); % 速度更新 vel(p,:) w * vel(p,:) ... c1 * r1 .* (pbest(p,:) - pos(p,:)) ... c2 * r2 .* (gbest - pos(p,:)); % 速度限幅 vel(p,:) max(min(vel(p,:), vmax), -vmax); % 位置更新 pos(p,:) pos(p,:) vel(p,:); % 位置边界修正 pos(p,:) max(min(pos(p,:), numCand), 1); end % 评价所有粒子 for p 1:popSize sel discreteSelect(pos(p,:), numCand, numSelect); cost calcCost(sel, demandCoord, demandWeight, candCoord, fixedCost); if cost pbestCost(p) pbestCost(p) cost; pbest(p,:) pos(p,:); end end % 更新全局最优 [minCost, minIdx] min(pbestCost); if minCost gbestCost gbestCost minCost; gbest pbest(minIdx, :); end bestCostHistory(iter) gbestCost; end % 输出最优方案 selFinal discreteSelect(gbest, numCand, numSelect); fprintf(最优中心编号: %s\n, mat2str(selFinal)); fprintf(最优总成本: %.2f\n, gbestCost);速度限幅这个动作别小看。如果不限幅粒子速度可能越飞越大位置在候选点编号之间剧烈震荡算法直接发散。我一般把最大速度设为位置取值范围宽度的0.5倍也就是0.5 * (numCand - 1)实际效果比较稳。3.5 结果可视化收敛曲线与选址地图迭代完之后最好把结果画出来。一方面是验证算法确实在收敛另一方面是给业务方看选址方案时一张图比一堆数字有说服力得多。绘图代码如下% 绘制收敛曲线 figure(Position, [100, 100, 1200, 500]); subplot(1, 2, 1); plot(1:maxIter, bestCostHistory, LineWidth, 2); xlabel(迭代次数); ylabel(最优成本); title(粒子群收敛曲线); grid on; % 绘制需求点、候选点和最终选中点 subplot(1, 2, 2); hold on; plot(demandCoord(:,1), demandCoord(:,2), b*, MarkerSize, 8); plot(candCoord(:,1), candCoord(:,2), ks, MarkerSize, 10, LineWidth, 2); plot(candCoord(selFinal,1), candCoord(selFinal,2), ro, MarkerSize, 14, LineWidth, 2); % 画需求点到最近中心的连线 for i 1:numDemand distVector sqrt((demandCoord(i,1) - candCoord(selFinal,1)).^2 ... (demandCoord(i,2) - candCoord(selFinal,2)).^2); [~, nearestIdx] min(distVector); nearestPoint candCoord(selFinal(nearestIdx), :); plot([demandCoord(i,1), nearestPoint(1)], ... [demandCoord(i,2), nearestPoint(2)], -, Color, [0.6, 0.6, 0.6]); end xlabel(X坐标); ylabel(Y坐标); legend(需求点, 候选点, 选中中心, 分配关系, Location, best); title(配送中心选址结果); grid on;运行之后你会看到一条逐渐下降并趋于稳定的收敛曲线以及一张包含需求点、候选点、选中点和连线的选址图。连线代表每个需求点的最近分配关系能直观看出哪些需求点被划给了哪个中心。4. 代码跑通之后的“三座山”参数、离散化、早熟4.1 惯性权重递减不是玄学是为了先粗后细很多粒子群教程直接把惯性权重设为固定值比如0.6。运行100次可能效果还行但一旦问题规模变了固定的w很容易要么前期收敛太慢要么后期震荡太凶。我习惯用线性递减w从0.9降到0.4。为什么算法前期需要较大的惯性权重让粒子保持较快的飞行速度在解空间里大范围探索避免过早扎进某个局部最优。到了后期粒子应该慢下来在当前最优解附近精细搜索所以w要小。线性递减是最简单的方案也是性价比最高的方案。你也可以尝试非线性递减比如按指数衰减但在我这个算例里线性递减和指数递减的结果差异非常小没必要为了花哨增加调参复杂度。如果你发现算法前期就卡住了可以把wStart提高到1.0如果后期震荡严重把wEnd降到0.2。这两个值得根据问题规模微调不要照搬所有参数。4.2 连续变量映射到离散选址最容易踩的重复坑前面提到过连续粒子维度取整后可能重叠导致实际选中的中心数不足。这个问题的根源在于PSO搜索空间是连续的而选址决策空间是离散的两者之间的映射不可能做到一一对应。我在discreteSelect里用“去重补齐”的方式处理。这种方法简单但也有一个副作用破坏了粒子原本携带的搜索方向信息。比如粒子本来想往编号2移动结果被随机补成了编号5那这次更新就相当于白跑了。为了缓解这个问题可以考虑在适应度函数中增加惩罚项强制重复编号产生很大的成本引导粒子避免重复。但罚函数法会大幅增加函数评估的复杂性小算例里没必要。更稳妥的做法是直接换二进制PSO粒子每一位对应一个候选点速度先映射成0或1的选择概率。不过二进制PSO对“选k个”的约束处理也比较麻烦通常要在解码时修复。说到底离散PSO本身就是工程化的算法没有完美方案能解决问题就行。4.3 早熟收敛的识别与补救粒子群最常见的失败模式是早熟收敛。现象是收敛曲线早早变得平直全局最优在几十次迭代内就不再下降但你知道这不是最优解。原因是种群多样性不足所有粒子都被吸引到当前全局最优附近失去了探索新区域的能力。识别早熟有几个办法一是看全局最优是否连续二三十次迭代没有任何变化二是看粒子位置之间的距离是否普遍小于某个阈值。如果粒子都聚在一起说明多样性已经耗尽。补救措施有三类按性价比排序对全局最优施加扰动每隔一定迭代次数在全局最优位置附近加入随机偏移再重新评价。这个实现最简单。重新初始化部分粒子将粒子群中适应度最差的一部分重新随机生成让它们重新出发探索。动态调整惯性权重检测到早熟时把w突然调高让粒子重新飞散。我在实际项目里用得最多的是“在全局最优上做扰动”每迭代20次把gbest按10%的幅度随机抖动一次。这样做实现简单也不太破坏已有的收敛趋势。5. 从仿真结果到真实决策谨慎对待“最优解”5.1 距离口径欧氏距离计算的成本误差有多大目前算例里的距离都是欧氏距离也就是直线距离。但真实物流中货物从仓库到客户走的是道路距离至少是直线距离的1.2到1.5倍如果遇到河流、山地、单行道绕行距离会更夸张。不同方向的道路密度不一样直线距离算出的成本比例和真实路况差距很大。如果只是做方案比选用欧氏距离误差尚可接受因为所有候选点都在同一套坐标系下相对优劣不会完全反转。但如果要把成本数字直接用来做投资估算就必须用路网距离。做法是把需求点和候选点导入地图工具取路网路径或者使用GIS软件的路网分析模块计算真实行驶距离。MATLAB里可以用Mapping Toolbox的distance函数做球面距离但还不是路网距离。更实际的做法是先跑一遍仿真确定候选方案集合再针对候选方案用路网距离做二次校验。5.2 需求量和成本数据的敏感性最优解是脆弱的粒子群给出的“最优解”是在特定数据和参数下成立的。真实业务中需求量经常变化新客户加入、老客户流失、大促前需求暴增这些都会让最优选址方案发生漂移。如果方案A比方案B只省了1%的成本但A对未来需求波动的适应性很差那实际决策更应该选B。我的习惯是跑完标准算例之后做一组敏感性分析分别把需求量上下浮动20%把固定成本上下浮动30%重复运行粒子群观察最优选中点是否变化。如果某些候选点在多次场景中反复出现说明这些点是稳健的可以优先考虑。如果每次结果都变说明系统对参数太敏感需要采集更准确的数据后再做决策。5.3 当问题规模变大粒子群还能用吗本文的算例只有30个需求点、6个候选点属于教学级别。真实场景可能是几百个需求点、几十个候选点。这时粒子群的粒子维度等于要选的k不算太高但每次适应度计算都要遍历所有需求点和所有候选点计算量会明显上升。好在MATLAB的向量化运算能扛住只要写好函数几百个需求点也只是一两秒的事。如果候选点数量非常大比如全省几百个乡镇街道作为候选点要用二进制编码的PSO或者改用遗传算法加局部搜索的组合策略。我的个人经验是先用K-means或贪心算法初选一批候选点缩小范围再用粒子群在缩小后的候选池里精细选点。这样既能控制计算量又比纯启发式结果更稳。我在实际项目里还习惯多跑几个随机种子比如分别用rng(1)到rng(10)各跑一遍记录每次选出的最优集合和成本。这样能看出算法稳定性。如果10次里8次选出的中心点相同只有成本略有波动那这个结果就可以拿去汇报。如果每次选点都不一样说明数据敏感性太高先回头查距离计算和参数设置而不是急着用结果。