经济的全球化和竞争的白热化迫使制造企业采用分布式布局，将工厂建设在离客户更近的地方，旨在实现快速响应市场需求、降低生产成本并缩短制造周期的目标[1]．但在此布局下，不同工厂发往同一客户的产品配送时长差异显著．现有生产计划编制中少有考虑产品配送过程，可能因配送能力受限、配送时间过长等原因，无法及时完成产品交付任务，致使生产计划失效、客户满意度降低、企业竞争力降低．与此相反，若能高效整合生产和配送过程，则最多降低20%的总运营成本[2]．因此，有必要将分布式生产中的产品配送过程纳入生产计划编制流程，开展分布式生产与配送的集成调度工作．针对相关问题，目前学术界有不少探讨．就研究问题而言，生产环境主要是分布式同构流水车间[3]或作业车间[4]，配送决策为路径规划或直达批量配送，还可能考虑共享配送资源[5]、不确定配送时长[6]、配送仓库[7]等其他特征，而目标函数基本与时间、成本或客户满意度相关．实际上，呈分布式布局的各车间的生产环境、生产资质大不相同；且在设备性能和产品工艺路线限制下，不同车间的生产运行状况也不同，故须围绕更真实的分布式生产场景，研究生产与配送的集成决策．从求解方法来看，鉴于问题的强NP难特征，目前主要采用启发式和元启发式算法进行大规模问题求解．当各车间生产环境异构、生产运行方式存在差异时，采用此类算法很难克服问题本身的复杂性，在编码和解码过程中也难以确保求解结果的近优性．相对而言，仿真优化不仅能模拟不同车间的动态行为属性，还能通过多次仿真运行进行算法迭代寻优，故而采用仿真优化来求解分布式生产与配送集成调度问题．因此，受到某易腐产品企业生产启发，聚焦于一类具有现实意义的分布式异构流水车间生产与配送集成调度问题，该问题允许部分阶段跨越．在此基础上，构建了融合多智能体系统(multi-agent system，MAS)和多目标头脑风暴优化(multi-objective brain storm optimization，MOBSO)算法的智能仿真优化框架，以探索分布式工厂生产与配送的同步决策方法．通过最晚送达时间和总成本的同步最小化，提升作业效率、降低作业成本．其中，分布式生产与配送过程的动态行为属性模拟和考虑问题特征的算法组件改进是关键．1 问题描述与数学建模1.1　问题描述呈分布式布局的多家工厂将联合完成一批生产订单．在生产过程中，部分订单无须在某些生产阶段进行加工，即存在工序跳跃的现象[8]．每个工厂按照流水车间组织生产，其中个性化定制产品在生产过程中部分阶段允许跨越．产品制造完成后由各工厂安排车辆，向客户进行产品配送．问题详细描述见图1，目标是最小化制造总成本及最晚送达时间．该问题具有以下现实性特征：a．各工厂异构，具有不同数量的生产阶段，且各阶段的加工工艺、生产能力均不同；b．各工厂呈流水车间布局，在产品制造过程中允许跳过部分生产阶段，不进行生产；c．每一订单仅由一家具有资质的工厂进行生产．10.13245/j.hust.250609.F001图1分布式生产与配送集成调度示意图该问题的决策变量有四个：订单在工厂间的分配、工厂内订单的加工顺序、订单的配送车辆安排及车辆针对多订单的配送路径规划．前两个决策是解决生产阶段的分布式异构流水车间调度问题，后两个决策是解决配送阶段的多车场多车辆路径规划问题．问题假设：a．工厂中所有工作中心在零时刻可用；b．每个工作中心同一时刻只能执行一个订单；c．每个订单同一时刻只能在一个工作中心内加工；d．每个订单一旦开始加工不允许中断；e．每个订单须由一辆车配送；f．每趟配送行程结束后，车辆须返回所属工厂；g．每辆车有相同的最大载重量限制；h．单个订单物料质量不能超过车辆最大载重量；i．每个订单至少有一家工厂具有生产资质．1.2　数学模型目标函数如下minTmax；(1)    min∑j∈J/{0}∑g∈Fj∑i∈Mg∑t∈TXjgitpjgiαg+βg∑g∈F∑h∈Kg∑j∈J/{0}∑k∈J/{0}Zghjkτjk+∑g∈F∑h∈Kg∑j∈J/{0}Zghj0σjg+∑g∈F∑h∈Kg∑k∈J/{0}Zgh0kσkg+∑j∈J/{0}CIj+∑j∈J/{0}CUj, (2)式中：Tmax为所有订单的最晚送达时间，作为优化目标，式(2)中总成本按照顺序包含了生产成本、配送成本、库存成本、惩罚成本四项，也作为优化目标；J为订单索引集合，J={0，1，…，j，k，…，n}，n为订单数量，0为虚订单索引；F为工厂索引集合，F={1，2，…，g，…，f}，f为工厂数量，Fj为具有加工订单j的工厂索引集合，Fj⊂F；M为所有工作中心类型的索引集合，M={1，2，…，i，m，…，|M|}；Mg为工厂g的工作中心类型的索引集合，Mg⊂M，且M1∪M2∪…∪Mg∪…∪Mf=M；K为全部工厂车辆的索引集合，K={1，2，…，h，…，v}，v为车辆总数；Kg为工厂g的车辆集合，Kg⊂K，且K1∪K2∪…∪Kg∪…∪Kf=K；T为工作中心订单加工顺序索引集合，T={1，2，…，t，…，n}；pjgi为订单j在工厂g工作中心i的加工时长；τjk为从订单j到订单k的配送时长；σjg为从工厂g到订单j的配送时长；αg，βg，γg分别为工厂g单位时间的生产成本、配送成本及库存成本；Xjgit为决策变量，若订单j在工厂g工作中心i中第t个加工则为1，否则为0；Zghjk为决策变量，若工厂g车辆h送完客户j后直送客户k则为1，否则为0；CIj和CUj分别为订单j的库存成本和惩罚成本．以下公式定义了两个目标的计算．Tmax≥Ahj    (∀h∈K,  ∀j∈J);Cij≥γg(Dh-Cji)-G(1-Ughj)(∀j∈J/{0},  g∈Fj,  h∈Kg,  i∈Mg);CUj≥ηωj(Ahj-bj)-G(1-Ughj)(∀j∈J/{0},  g∈Fj,  h∈Kg);CUj≥ηωj(aj-Ahj)-G(1-Ughj)(∀j∈J/{0},  g∈Fj,  h∈Kg),式中：Ahj为车辆h到达客户j的时间；Cij为订单j在工作中心类型i上完成加工时间；Dh为车辆h离开工厂的时间；G为一个足够大的正数；η为提前期和延迟期的单位时间惩罚成本；aj，bj，ωj分别为订单j的最早交货期、最晚交货期及权重；Ughj为决策变量，若订单j由工厂g的车辆h进行配送则为1，否则为0．订单分配约束：除虚订单外，每个订单须由一家工厂进行加工，并且在该工厂的对应工作中心上加工一次．∑g∈FjYjg=1    (∀j∈J/0);∑t∈TiXjgit=Yjg(∀j∈J/{0},  g∈Fj,  s∈Oj,  i∈Mg,ljsi=1),式中：Yjg为决策变量，若订单j分配给工厂g加工则为1，否则为0；Oj为订单j的工序索引集合，索引为s；ljsi为0-1参数，若订单j的第s道工序需要类型为i的工作中心则为1，否则为0．加工顺序约束：一个订单的一个工序必须分配到一个工厂的对应工作中心的队列中，并且按照订单工序顺序进行加工．∑g∈Fj∑t∈TiXjgit=1    (∀j∈J/{0},i∈M);∑j∈J/{0}∑s∈Oj,ljsi=1Xjgit≤1(∀g∈Fj,  i∈Mg,  t∈T);∑j∈J/{0}∑s∈Oj,ljsi=1Xjgi(t-1)≤∑j∈J/{0}∑s∈Oi,ljsi=1Xjgit(∀g∈Fj,  i∈Mg,  t∈T).工序定时约束：订单每道工序的开始时间晚于前一工序结束时间，其加工时长严格符合工艺要求．其中，Bgit和Fgit分别为在工厂g中工作中心i的第t个加工任务的开始时间和结束时间．Bgi(t+1)≤Fgit    (∀g∈F,  i∈Mg,  t∈T,t≠T);Fgit=Bgit+∑j∈J/{0}Xjgitpjgi    (∀g∈Fj,  i∈Mg,t∈T);Cji≥Fgit-G(1-Xjgit)(∀j∈J/{0},  g∈Fj,  i∈Mg,  t∈T);Cji≤Fgit-G(1-Xjgit)(∀j∈J/{0},  g∈Fj,  i∈Mg,  t∈T);Bgmt≥Cji-G(1-Xjgmt)(∀j∈J/{0},  g∈Fj,  t∈T,  i,  m∈Mg).车辆数量约束：分配给工厂的每个订单都有车辆负责承运．∑h∈KgUghj=Yjg    (∀j∈J/{0},  g∈Fj).配送路径约束：为完成所分配的订单配送任务，每辆车从所属工厂出发、依序配送、返回所属工厂．∑k∈JZghjk=Ughj(∀j∈J/{0},  j≠k,  g∈Fj⋂Fk,  h∈Kg);∑k∈JZghkj=Ughj(∀j∈J/{0},  j≠k,  g∈Fj⋂Fk,  h∈Kg);∑k∈JZghjk=Ughj(∀j∈J/{0},  j≠k,  g∈Fj⋂Fk,  h∈Kg);∑k∈JZghkj=Ughj(∀j∈J/{0},  j≠k,  g∈Fj⋂Fk,  h∈Kg);Zghjk+Zghkj≤Yjg(∀j,  k∈J/{0},  g∈Fj⋂F,  h∈Kg);∑h∈KgZghj0≤1    (∀j∈J/{0},  g∈Fj);∑h∈KgZgh0j≤1    (∀j∈J/{0},  g∈Fj).车辆载重约束：每辆车所配送订单的物料总质量不超过车辆的最大载重量，其中，dj和dmax分别为订单j的物料质量及每台配送车辆的最大载重量，对于任意订单j，dj≤dmax．          ∑h∈KgZghjkdj≤dmax(∀j, k∈J/{0},  k≠j, g∈Fj⋂Fk).配送定时约束：每一订单所在车辆开始配送的时间须晚于订单生产完成时间；配送过程中须在前一订单配送结束后，再开始后续订单配送．Dh≥Cji-G+G∑k∈JZghkj(∀j∈J/{0},  g∈Fj,  h∈Kg,  i∈Mg);Ahj≥Dh+σjg-G(1-Zgh0j)(∀j∈J/{0},  g∈Fj,  h∈Kg);Ahk≥Ahj+τjk-G(1-Zghjk)(∀j, k∈J/{0},  j≠k,  g∈Fj⋂Fk,  h∈Kg).2 基于智能仿真优化的头脑风暴算法上述数学模型用GAMS软件对小规模问题采用7个订单、3个工厂的EP1进行求解．根据求解器得到的一个非支配解，在仿真优化模型中固定所有的决策变量，结果表明生成的调度方案完全一致，验证了数学模型的准确性．由于问题是NP难问题，无法在合理的时间内求解中大规模问题，因此采用优化算法进行求解．由于本问题中多工厂的布局异构、资质各异、生产过程中跨越阶段各不相同，因此以仿真软件为开发平台，根据建立的数学模型构建多智能体系统，对异构车间的动态行为属性进行高保真模拟．最新的此类研究中利用头脑风暴算法进行优化，取得了出色的效果[9]，故提出融合多智能体系统与多目标头脑风暴优化的仿真优化框架(见图2)，通过强强结合，发挥各自优势．多目标头脑风暴优化的主要改进点包括成本目标导向的初始化启发式规则、基于编码差异性与支配关系的集群划分方法、基于关键工厂的邻域搜索算子，同时设计了基于概率的集群中心替换策略和基于问题特征的交叉/变异与修复机制．10.13245/j.hust.250609.F002图2基于智能仿真优化的头脑风暴法框架2.1　基于多智能体系统的动态行为属性模拟由于各工厂布局异构、资质不同、运行方式不同，因此采用三个启发式规则，限制不同工厂的运行流程、车辆分配和车辆配送属性，模拟异构环境下车间运行的动态行为属性，多智能体系统的时序逻辑见图3．10.13245/j.hust.250609.F003图3基于多规则的仿真时序逻辑a．工作中心排序规则：基于先到先服务，进行各工作中心工序分配．值得注意的是，编码中仅限定了订单在工厂的投产顺序，而订单在不同工作中心的加工顺序则根据该规则确定．b．车辆分配规则：基于订单加工完成时间，安排最早空闲车辆配送，若当前车辆加入该订单后达到最大载重，则该订单分配到下一辆车．c．路径规划规则：基于最近邻法，确定下一个配送的订单，直至所分配订单完成为止．由此可见：当订单工厂分配及厂内投产顺序已知时，依据工作中心排序规则，松弛待执行工序操作与具体订单工艺路线的关系，能够驱动异构工厂正常运行．在工厂内订单加工完成时间未知的情况下，利用车辆分配规则可优先将加工完成时间相近的订单组批，实现车辆的合理配送；而当车辆批次订单数不确定时，若对配送路径进行编码，会导致空间增大．此时，利用基于车辆动态位置的路径规划规则进行配送路径决策，可保证车辆有序配送．2.2　成本目标导向的初始化启发式规则针对该问题，采用向量组编码方式，一个染色体就是一个方案解．如{{5，2，1}，{3，6}，{4，7，8}}表示8个订单分到3家工厂及在各工厂的投产顺序，可知第一家工厂共执行3个订单，先订单5再订单2和1．改善种群初始化性能，有利于提高种群多样性、加快算法收敛．因此，采用混合初始化策略，用随机方法产生一半初始解，由启发式规则产生另一半初始解．混合初始化策略的有效性在实验部分进行了验证．考虑到每个订单在各工厂的制造成本及从各工厂送达客户的配送成本均不同，由此计算每个订单在各工厂的生产/配送的估计成本，将该订单优先分给估计成本最低的工厂，如下式所示，g=argming∈Fαg∑s∈Oj,i∈Mg,ljsi=1pjgi+2βgσjg(∀j∈J), (3)括号中第一项是订单在工厂的生产总成本，第二项是工厂和订单之间的往返配送成本．但因为订单配送顺序未知，所以订单之间的配送在启发式规则中不考虑．2.3　种群的集群划分及新种群生成策略为实现相似解的聚类，并确定各个集群的中心，提出了一种基于编码差异性与支配关系的集群划分方法．由于待求解问题具有多模态特性，即存在适应度相同但解不同的情况，为量化表征解之间的差异，因此提出一种基于染色体编码差异性的K均值聚类方法，用于将相似解聚集到同一集群．定义Yjg1和Yjg2为两个解中工厂分配的决策变量，统计订单的工厂分配决策差异，即度量解的差异化程度ζ，ζ=∑g∈F∑j∈J/{0}|Yjg1-Yjg2|. (4)聚类的具体过程如下：a．须要指定集群聚类的集群数量K，并随机选择K个对象作为初始的聚类中心；b．将数据集中每个对象分配给与其差异化最小的聚类中心，形成K个聚类组；c．重新计算每个类中的聚类中心；d．重复步骤b和c，直至聚类中心不再变化．此外，经典头脑风暴算法在单目标优化中，集群中心只有一个，即适应度最优的解；扩展到多目标优化，集群中心由该集群中所有的非支配解构成(见图4)．10.13245/j.hust.250609.F004图4基于编码差异性与支配关系的集群划分通过种群的集群划分，集群内部的解具有较大相似性，而集群之间的解差异相对较大．新种群生成策略考虑编码差异性，具体步骤如下．a．产生随机数r1∈[0，1]，若r1小于单集群概率Pb，则执行步骤b，否则跳转到步骤c．b．随机选择一个集群S，产生随机数r2∈[0，1]，若r2小于选择中心概率Pc，则采用轮盘赌从S的集群中心选择一个解π，否则采用轮盘赌从S中选择一个解π．对π进行变异或邻域搜索操作，产生新解，并添加入新种群．c．随机选择两个集群S1和S2，产生随机数r3∈[0，1]，若r3小于选择中心概率Pc，则采用轮盘赌从S1和S2的集群中心各选择一个解，否则采用轮盘赌从S1和S2中各选择一个解，对选出的两个解进行交叉操作，产生两个新解．d．将产生的新解全部加入新种群，若新种群不小于规定种群大小，则结束种群生成，否则跳转到步骤a．2.4　交叉、变异与修复经典头脑风暴算法适用于连续优化问题，而本研究问题属于离散优化问题，其原始算子无法直接应用；因此，针对性地采用OX/PBX交叉算子、单点/两点变异算子及基于关键工厂的邻域搜索算子．其中OX/PBX交叉算子是随机选择离散几个或一段连续的基因位，将父代1的相应位置基因直接继承到子代中，将父代2的余下基因顺次填充到子代中．单点变异算子是随机选择一个并插入到其他位置中，而两点变异算子是随机选择两个基因进行互换．基于关键工厂的邻域搜索算子的具体步骤如下：a．获取当前解的最早、最晚完成配送任务的工厂索引fmax和fmin；b．找到fmax对应工厂所有订单，并从中找到在fmin对应工厂中生产与配送成本之和最小的订单，将该订单插入fmin对应的订单序列的随机位置．鉴于交叉、变异和邻域搜索算子未考虑工厂生产资质的差异，可能产生不可行解，为此设计了一种染色体修复算法，具体步骤如下：a．检查每个工厂的订单序列中是否存在无法加工的订单；b．若存在无法加工订单，则随机选择一个能加工的工厂，并将该订单插入该工厂对应订单序列的随机位置．由于问题的可行解空间极大，为保证种群的稳定性和多样性，因此选择支配层级高或同一层级但拥挤距离大的解进入下一代种群．3 实验与分析为验证所提算法的有效性和优越性，设计了算法参数校验、改进算子性能验证和算法对比三类实验．所有算法均采用Java语言编写，在配备Core(TM) i5-11400H CPU和16 GiB RAM的计算机上执行．算法性能评价指标采用超体积率(H)、逆世代距离(I)和世代距离(G)，算法终止时间根据订单数量n进行设置，当n＞20时设置为300 s，反之设置为180 s．3.1　算例生成参考实际生产实例，获取150个不同订单，这些订单涵盖6种典型工艺路线，据此构建30个测试案例，分别记为EP1～EP30．案例的参数设置如下：订单权重、订单数量、订单工序数、工厂数量、加工时间(h)、订单质量(t)、横纵坐标(km)、最早交货期(h)、最晚交货期(h)分别在U[0，1]，U[7，150]，U[3，6]，U[3，7]，U[Ts-2，Ts+2]，U[1，30]，U[1，30]，U[200，800]，U[60，Tsum]-U[10，70]/2，U[60，Tsum]+U[70，130]/2取值．其中：Ts为各订单的标准工序时间；Tsum为订单所有工序的标准加工时间总和；U[x，y]为一个在区间[x，y]上均匀分布的随机变量．3.2　参数校验为了评估所提出的多目标头脑风暴优化算法的有效性和优越性，将其与用于解决同类问题的其他算法进行性能对比，对比算法包括第二代非支配遗传算法(NSGA-Ⅱ)[10]、基于分解的多目标进化算法(MOEA/D)[11]、多目标教与学算法(MTLBO)[5]．此外，由于问题特征及编码方式不同，对比算法进行了相应修改．实验选取包含120个订单、7个工厂的测试案例EP30，该案例规模较大且具有典型特征．对其进行参数调优，有助于评估算法在实际应用中的性能表现．H是评价收敛性和多样性的综合性指标，无须知道真实前沿或参考集，采用H作为响应．根据田口法设计16种参数组合，每种组合重复运行20次后取H均值．对比算法均采用参数正交试验进行参数校验，并通过算法的参数主效应图选择H均值最大的参数水平．多目标头脑风暴优化算法的参数主效应分析见图5，所有对比算法的参数设置如下．多目标头脑风暴优化算法：种群大小Ps=60，集群数量K=10，随机替换集群中心的概率Pa=0.6，Pb=0.4，Pc=0.4．MOEA/D算法：子问题数为40，邻居数量为4．NSGA-Ⅱ算法：种群大小为20，交叉概率为0.4，变异概率为0.2．MTLBO算法的种群大小为60，教学比例为0.4．10.13245/j.hust.250609.F005图5多目标头脑风暴优化算法参数主效应图3.3　改进算子有效性校验针对本研究的问题，在所提算法中进行了三方面改进，分别为基于编码差异性与支配关系进行集群划分、成本目标导向的初始化启发式规则和基于关键工厂的邻域搜索算子．为检验所提出的基于编码差异性的聚类(记为KC)的有效性，将其与基于支配关系的聚类(记为DC)[9]、基于目标函数值的聚类(记为SC)、基于随机的聚类(记为RC)相对比．聚类方法校验结果见图6，可见所提出的聚类方法的鲁棒性更高，基于编码差异性的聚类能够将差异性更低、相似性更高的解置入同一集群，通过相似解的共同演化和差异解间的基因交流，可以有效促进算法的局部开发和全局探索．10.13245/j.hust.250609.F006图6聚类方法校验将多目标头脑风暴优化算法中移除成本目标导向的初始化启发式规则(HI)后的算法表示为MOBSO_HI，将多目标头脑风暴优化算法中移除基于关键工厂的邻域搜索算子(MS)后的算法表示为MOBSO_MS．实验结果如图7所示，在30个测试案例中，多目标头脑风暴优化算法在三个指标的均值与方差方面都显著优于其他三个算法，分析可知：初始化启发式规则明显提高了算法初始种群的质量，加快了算法收敛速度，而邻域搜索算子增强了算法的局部搜索能力．10.13245/j.hust.250609.F007图7算法组件校验3.4　算法有效性验证根据四种对比算法在30个测试案例中得到的评价指标，绘制箱线图见图8．10.13245/j.hust.250609.F008图8算法评价指标对比采用Kruskal-Wallis检验方法结果显示：在H，G，I三个指标上，P值都小于0.001，说明四种算法在性能上差异有统计学意义．在指标H方面，多目标头脑风暴优化算法仅在EP20案例中未显著优于NSGA-II，MOEA/D和MOTLBO算法，在其余案例中均表现更优；在指标G方面，除EP16和EP27外，多目标头脑风暴优化算法优势明显；在指标I方面，除EP4，EP11，EP12和EP27外，多目标头脑风暴优化算法更具优势；而在方差表现上，多目标头脑风暴优化算法始终最佳，稳定性远超其他算法．综上所述，实验结果表明多目标头脑风暴优化算法相较于其他算法，在性能上更为稳定有效．深入分析算法机理可知，成本目标导向的初始化启发式规则提升了初始种群质量，基于编码差异性与支配关系的集群划分方法平衡了算法的局部开发与全局探索能力，而基于关键工厂的邻域搜索算子则进一步强化了算法的局部开发能力．4 结语针对分布式异构、带工序跳越的流水车间生产与配送集成调度问题，以最晚配送完成时间、总成本最小化为目标，建立混合整数线性规划模型，并提出基于多智能体系统和多目标头脑风暴算法的智能仿真优化框架求解．利用多智能体系统模拟生产运输过程中的动态行为属性，根据编码差异性与支配关系进行集群划分，利用成本最小化目标建立初始化启发式规则，基于关键工厂设计邻域搜索算子，实现算法局部开发和全局探索的平衡，以实现高效、快速求解问题．未来研究将考虑生产和配送中的不确定性因素，如订单临时到达、机器故障等，研究分布式生产与配送的鲁棒集成调度问题．