原材料订购多目标优化模型及求解算法

左 霞 封汉颍

（南通理工学院 基础教学学院，江苏 南通 226000）

1 引言

随着我国市场经济的持续发展，原材料的订购成本控制越来越被生产企业关注，合理的订购不仅能够降低企业经营风险，还能提高企业竞争能力。实际采购中，订购方案除了追求最经济的目标外，还需要考虑对多种原材料的选择、供货商的数量等等，因此优化一个多目标的原材料订购方案尤为必要。

多目标优化在工程应用、推荐系统、物流配送、路径规划等各个领域中都普遍存在，但是多个性能目标之间常存在不可比较性和相互矛盾的现象，不一定在所有目标上都有最优解，即所有目标函数必须达到最优。文献[1]用遗传算法得到Pareto最优解后用TOPSIS方法生成最优方案；文献[2]求出Pareto最优解集供决策者和设计者选择；文献[3]对多目标优化模型进行求解，得到了Pareto前沿解集，提出了基于加权相对距离的决策方法，得到了最优生产方案；文献[4]给出了对已有的Pareto最优解用R方法评价，求得最佳折中方案。

解决多目标优化的方法有MOEA算法、NSGA-II算法、粒子群算法，以及其他智能算法如遗传算法、蚁群算法、人工蜂群算法等。其中普通的粒子群算法的粒子初始位置、更新速度都是连续函数，与之对应，位置和速度更新均为离散值的算法是离散PSO算法。BPSO是在离散粒子群算法基础上，约定位置向量、速度向量均由0、1值构成；BPSO有很强全局搜索能力，但不能收敛于全局最优值，且随着算法迭代搜索随机性越来越强，缺乏后期的局部搜索能力。NSGA-II算法的特点是Pareto最优解分布更均匀，因此，本文采用这两种算法对实例进行仿真实验，并对结果进行对比。

2 问题描述及其数学模型的构建

2021年高教社杯数学建模竞赛C题，是以某生产企业订购与运输原材料为背景的竞赛题。其中制定的订购方案，有要求供应商最少、有要求最经济，以及对于多种原材料有选择性地订购等。在实际订购中，决策者往往希望这三个要求都能满足，故本文把这三问稍加改变和聚合，得到如下问题，其表述为：某生产企业所用原材料总体可分为A，B，C三种类型，该企业每年按48周安排生产，需要提前制定24周的原材料订购，即根据产能要求确定需要订购的原材料供应商（称为“供应商”）和相应每周的原材料订购数量（称为“订货量”）。该企业每周的产能为2.82万立方米，每立方米产品需消耗A类原材料0.6立方米，或B类原材料0.66立方米，或C类原材料0.72立方米。由于原材料的特殊性，供应商不能保证严格按订货量供货，实际供货量可能多于或少于订货量。为了保证正常生产的需要，该企业要尽可能保持不少于满足两周生产需求的原材料库存量，为此该企业对供应商实际提供的原材料总是全部收购。原材料的采购费用对企业的生产收益有直接的影响，实际中A类和B类原材料的采购单价分别比C类原材料高20%和10%。三类原材料运输和储存的单位费用相同。现在有50家最重要的供应商，决策者希望选择最少的供应商，根据前5年的供货量制定未来24周每周最经济的原材料订购方案，同时该企业为了控制生产成本，计划尽量少地采购A类和尽量多地采购C类原材料。

2.1 建立目标函数的模型

根据决策者订购的要求，按照不同的优化目标建立相应的模型，具体分为以下三步：

第一步，把选择最少的供应商作为目标函数。设dij表示第i个供货商在第j周的供货量，由dij构成了一个供货量的矩阵D。设决策变量为xij，当xij=1，代表第i个供货商在第j周被选上，当xij=0，代表第i个供货商在第j周没有被选上，由xij构成决策变量矩阵X，现在建立如下数学表达式：

其中T为所有供货商对应的指标集，N表示计划订购的总周次，sgn(x)为符号函数。

第二步，把原材料订购费用作为目标函数，求最经济的订购方案。设mij表示第i个供货商在第j周的订货量，因为该企业对供应商实际提供的原材料总是全部收购，如果第i个供货商在第j周被选上，则mij=dij，如果第i个供货商在第j周没有被选上，则mij=0。根据前面的决策变量xij的定义，有mij=xij×dij。因为三类原材料运输和储存的单位费用相同，不妨设运输单位费用为a,储存的单位费用为b，假设a+b与原材料C的单价比是1:1，若原材料C的单价为1，原材料A和原材料B的单价分别为1.2和1.1，则M家供货商在N周内原材料订购的费用为：

其中TA，TB，TC是原材料A，B，C三类供货商集合分别对应的指标集。

第三步，以尽量少地采购A类和尽量多地采购C类为目标函数。通过增加权重转化成如下的单目标函数：

2.2 写出约束条件

因为总的供货量会在转运的过程中有损耗，根据近5年8家转运商的相关数据，平均的损耗率是1%，所以这里假设生产企业的接收量是供货量的99%。根据题目的说明，为了保障企业正常生产的需要，该企业要尽可能保持不少于两周生产需求的原材料库存量，企业每周产量是28200立方米，两周是56400立方米。按1%的损耗率得到约束条件为：

2.3 整体模型

记A=(aij)M×N，设B=[28484.8×2,28484.8×3,…,28484.8×25]T，B中第j个分量为Bj。再记A=(A1,A2,…,AN),X=(X1,X2,…,XN)，则约束条件表达式(2)整理为：

综上，该模型为非线性多目标0-1规划模型：minf=min(z1,z2,z3)，其中：

3 V-BPSO算法和NSGA-II算法

3.1 V型二进制粒子群算法（V-BPSO）

3.1.1 二进制粒子群算法（BPSO）

二进制粒子群算法（BPSO）最初是在1997年由J.Kennedy和R.C.Eberhart设计，是在离散粒子群算法基础上，约定位置向量、速度向量均由0、1值构成。算法流程如下：

初始化粒子位置：按一定策略，生成二进制编码；

第i个粒子的位置xi=(xi1,xi2,…,xiN)T，速度为vi=(vi1,vi2,…,viN)T。

速度更新公式：

位置的更新公式：

概率映射方式，采用sigmoid函数将速度映射到[0,1]区间作为概率，这个概率就是粒子下一步取值为1的概率；

位置变化的绝对概率：

当前位置为0变化为1，当前为1变化为0，这二者被称为绝对变化；概率表示为：

公式(6)中，i=1,2,…,M，M为该群体中的粒子总数，其中ω为非负的惯性因子，学习因子c1和c2是[0,2]之间的常数；rand1()和rand2()是介于[0,1]之间的随机数；Pbestid是粒子i在第d维的个体极值点的位置；Gbestd是整个种群在第d维的全局极值点的位置。

3.1.2 V型二进制粒子群算法（VBPSO）

BPSO有很强全局搜索能力，但不能收敛于全局最优值，且随着算法迭代搜索随机性越来越强，缺乏后期的局部搜索能力，因此有学者对BPSO算法进一步改进，Seyedali Mirjalili，Andrew Lewis提出二元粒子的S型和V型传递函数[5]，通过引入一个新的传递函数族来改进BPSO算法，在本文中，采用的位置更新函数公式为：

这里的(t)和(t)表示粒子i在第k维迭代t处的位置和速度，((t))-1是(t)的补码，这种方法的优点是这个传递函数不会强制粒子取0或1的值，它们鼓励粒子在速度值较低时停留在当前位置，或者在速度值较高时切换到它们的补码，（9）式被称为V型传递函数，粒子的位置更新图如图1所示。

图1 粒子的位置更新

3.2 非支配排序遗传算法（NSGA-II）

非支配排序遗传算法[6]是解决多目标优化的一种算法，它在选择、交叉、变异算子与普通的遗传算法相同，在这基础上，非支配排序遗传算法增加了两个部分，一个是非支配排序，另一个是定义拥挤距离及拥挤比较算子。

（1）非支配排序

在非支配排序算法中，每个解都会被分配两个参数np和Sp，其中np表示支配第p个解的个数，Sp为种群中被第p个解支配的个体集合。非支配排序算法的主要步骤为：

Step1：对于种群中np=0的个体，保存在当前集合F1中，并记其支配等级为1；

Step2：对于F1中的个体i，其所支配的个体集合为Si，遍历Si中的每个个体p，执行np=np-1，如果np=0则将个体i保存在集合H中，并记其支配等级为2；

Step3：以H作为当前集合，重复上述步骤，找到所有的非支配层。

（2）拥挤距离

拥挤距离描述群体个体拥挤程度（密度）和目标函数值有关。计算方法如下：

首先设个体的拥挤距离用d表示，个体i的拥挤距离用di表示，设置di=0，再设fm为目标函数，m=1,2,..M，将fm升序排列，第一位和最后一位的拥挤距离设为d1=d L=∞，对于非边界个体i拥挤距离度量值的计算公式为：

（3）拥挤比较算子

通过非支配排序以及拥挤度计算之后，种群中的每个个体n都得到两个属性：非支配序nrank和拥挤度nd。定义拥挤度比较算子为≥,个体优劣的比较依据为：i≥j，即个体i优于个体j，当且仅当irank≤jrank且id≥jd。

4 模型的求解及Pareto前沿图

本次求解实验基于MATLAB-R2020b进行编程，运行环境为：Windows 10(64位操作系统)，内存16GB处理器为Intel(R)Core(TM)i7-10510U CPU@1.80GHz 2.30 GHz。

首先设置VBPSO算法的参数，种群规模Np=150，最大迭代次数maxgen=500，惯性系数的初始值c1=c2=w=2，速度最大值vmax=6，最大惯性权重ωmax=0.9，最小惯性权重ωmin=0.4，其次对于约束条件利用外点法设置惩罚函数。实验数据见参考文献[7]，运行算法所得到Pareto前沿图如图2、3所示。

图2 VBPSO Prareto前沿

图3 NSGA-IIPrareto前沿

通过观察Pareto最优解可知，f2和f3两个优化目标之间存在着明显的约束关系，当f2减小时，f3是逐渐增大的，从而可证明本文所选目标函数符合多目标优化目标函数要求。这两种算法相比，VBPSO运行的速度比NSGA-II快，但是后者的Pareto前沿的点分布更均匀，计算更准确。

5 用R方法选择最佳解

在求得多目标优化的Pareto最优解后，决策者可以从可用的Pareto最优解方案中找到最佳替代解决方案。对于每一个Pareto最优解，优化目标之间的折中是不同的，这里用R方法寻找这种最佳的折中解决方案，具体的步骤如下：

步骤1：列一个决策表，包含与目标相对应的Pareto最优解方案的性能数据。

步骤2：根据决策者认为的1，2，3等的重要性对目标进行排名。如果两个或多个目标被认为同等重要，则为这些目标分配平均排名。

步骤3：根据目标相关的性能数据，按1，2，3等对Pareto最优解方案进行排名。如果两个或者多个方案具有相应于目标的相同值，则将平均排名分配给这些方案。

步骤4：将分配给目标等级和Pareto最优解方案转换为相应的权重。转换的公式如下：

wj为目标/备选的权重，rk为目标/备选的排序，n为目标/备选的数量。

步骤5：通过将目标权重与Pareto最优解方案的相应权重相加，计算各方案综合得分。

步骤6：根据综合得分，具有最高综合得分的方案被认为是最佳妥协Pareto最优解。

本模型由VBPSO求得Pareto最优解，共15个备选订购方案，对每个方案的三个目标函数值按数值的大小进行排序，因为f1数值相等，故只需要对f2和f3排序，同时按重要性排序f2＞f3＞f1，利用公式（12）转化成相应的权重。

通过目标权重与备选方案的相应权重相乘后再求和，得到备选方案的综合得分。其中方案4的综合得分最高，认为它是最佳折中的订购方案，综合得分见表1。

表1 各方案综合得分

6 结语

本文基于生产企业订购原材料的三个需求目标和该企业对供应商实际提供的原材料总是全部收购的特点，建立原材料订购方案的模型，最后转化成非线性的0-1规划问题。该模型采用V型的二进制粒子群算法（V-BPSO），引入一个新的函数更新粒子的位置，改进粒子群算法容易陷入局部最小的缺点，同时也提高收敛速度。用该算法对原材料订购模型求解，找到Pareto前沿，用R方法对备用方案评价，得到最佳折中解。因此该模型具有一定的理论价值和实践价值。