原材料訂購(gòu)多目標(biāo)優(yōu)化模型及求解算法

左 霞 封漢潁

（南通理工學(xué)院 基礎(chǔ)教學(xué)學(xué)院，江蘇 南通 226000）

1 引言

隨著我國(guó)市場(chǎng)經(jīng)濟(jì)的持續(xù)發(fā)展，原材料的訂購(gòu)成本控制越來(lái)越被生產(chǎn)企業(yè)關(guān)注，合理的訂購(gòu)不僅能夠降低企業(yè)經(jīng)營(yíng)風(fēng)險(xiǎn)，還能提高企業(yè)競(jìng)爭(zhēng)能力。實(shí)際采購(gòu)中，訂購(gòu)方案除了追求最經(jīng)濟(jì)的目標(biāo)外，還需要考慮對(duì)多種原材料的選擇、供貨商的數(shù)量等等，因此優(yōu)化一個(gè)多目標(biāo)的原材料訂購(gòu)方案尤為必要。

多目標(biāo)優(yōu)化在工程應(yīng)用、推薦系統(tǒng)、物流配送、路徑規(guī)劃等各個(gè)領(lǐng)域中都普遍存在，但是多個(gè)性能目標(biāo)之間常存在不可比較性和相互矛盾的現(xiàn)象，不一定在所有目標(biāo)上都有最優(yōu)解，即所有目標(biāo)函數(shù)必須達(dá)到最優(yōu)。文獻(xiàn)[1]用遺傳算法得到Pareto最優(yōu)解后用TOPSIS方法生成最優(yōu)方案；文獻(xiàn)[2]求出Pareto最優(yōu)解集供決策者和設(shè)計(jì)者選擇；文獻(xiàn)[3]對(duì)多目標(biāo)優(yōu)化模型進(jìn)行求解，得到了Pareto前沿解集，提出了基于加權(quán)相對(duì)距離的決策方法，得到了最優(yōu)生產(chǎn)方案；文獻(xiàn)[4]給出了對(duì)已有的Pareto最優(yōu)解用R方法評(píng)價(jià)，求得最佳折中方案。

解決多目標(biāo)優(yōu)化的方法有MOEA算法、NSGA-II算法、粒子群算法，以及其他智能算法如遺傳算法、蟻群算法、人工蜂群算法等。其中普通的粒子群算法的粒子初始位置、更新速度都是連續(xù)函數(shù)，與之對(duì)應(yīng)，位置和速度更新均為離散值的算法是離散PSO算法。BPSO是在離散粒子群算法基礎(chǔ)上，約定位置向量、速度向量均由0、1值構(gòu)成；BPSO有很強(qiáng)全局搜索能力，但不能收斂于全局最優(yōu)值，且隨著算法迭代搜索隨機(jī)性越來(lái)越強(qiáng)，缺乏后期的局部搜索能力。NSGA-II算法的特點(diǎn)是Pareto最優(yōu)解分布更均勻，因此，本文采用這兩種算法對(duì)實(shí)例進(jìn)行仿真實(shí)驗(yàn)，并對(duì)結(jié)果進(jìn)行對(duì)比。

2 問(wèn)題描述及其數(shù)學(xué)模型的構(gòu)建

2021年高教社杯數(shù)學(xué)建模競(jìng)賽C題，是以某生產(chǎn)企業(yè)訂購(gòu)與運(yùn)輸原材料為背景的競(jìng)賽題。其中制定的訂購(gòu)方案，有要求供應(yīng)商最少、有要求最經(jīng)濟(jì)，以及對(duì)于多種原材料有選擇性地訂購(gòu)等。在實(shí)際訂購(gòu)中，決策者往往希望這三個(gè)要求都能滿足，故本文把這三問(wèn)稍加改變和聚合，得到如下問(wèn)題，其表述為：某生產(chǎn)企業(yè)所用原材料總體可分為A，B，C三種類型，該企業(yè)每年按48周安排生產(chǎn)，需要提前制定24周的原材料訂購(gòu)，即根據(jù)產(chǎn)能要求確定需要訂購(gòu)的原材料供應(yīng)商（稱為“供應(yīng)商”）和相應(yīng)每周的原材料訂購(gòu)數(shù)量（稱為“訂貨量”）。該企業(yè)每周的產(chǎn)能為2.82萬(wàn)立方米，每立方米產(chǎn)品需消耗A類原材料0.6立方米，或B類原材料0.66立方米，或C類原材料0.72立方米。由于原材料的特殊性，供應(yīng)商不能保證嚴(yán)格按訂貨量供貨，實(shí)際供貨量可能多于或少于訂貨量。為了保證正常生產(chǎn)的需要，該企業(yè)要盡可能保持不少于滿足兩周生產(chǎn)需求的原材料庫(kù)存量，為此該企業(yè)對(duì)供應(yīng)商實(shí)際提供的原材料總是全部收購(gòu)。原材料的采購(gòu)費(fèi)用對(duì)企業(yè)的生產(chǎn)收益有直接的影響，實(shí)際中A類和B類原材料的采購(gòu)單價(jià)分別比C類原材料高20%和10%。三類原材料運(yùn)輸和儲(chǔ)存的單位費(fèi)用相同?，F(xiàn)在有50家最重要的供應(yīng)商，決策者希望選擇最少的供應(yīng)商，根據(jù)前5年的供貨量制定未來(lái)24周每周最經(jīng)濟(jì)的原材料訂購(gòu)方案，同時(shí)該企業(yè)為了控制生產(chǎn)成本，計(jì)劃盡量少地采購(gòu)A類和盡量多地采購(gòu)C類原材料。

2.1 建立目標(biāo)函數(shù)的模型

根據(jù)決策者訂購(gòu)的要求，按照不同的優(yōu)化目標(biāo)建立相應(yīng)的模型，具體分為以下三步：

第一步，把選擇最少的供應(yīng)商作為目標(biāo)函數(shù)。設(shè)dij表示第i個(gè)供貨商在第j周的供貨量，由dij構(gòu)成了一個(gè)供貨量的矩陣D。設(shè)決策變量為xij，當(dāng)xij=1，代表第i個(gè)供貨商在第j周被選上，當(dāng)xij=0，代表第i個(gè)供貨商在第j周沒(méi)有被選上，由xij構(gòu)成決策變量矩陣X，現(xiàn)在建立如下數(shù)學(xué)表達(dá)式：

其中T為所有供貨商對(duì)應(yīng)的指標(biāo)集，N表示計(jì)劃訂購(gòu)的總周次，sgn(x)為符號(hào)函數(shù)。

第二步，把原材料訂購(gòu)費(fèi)用作為目標(biāo)函數(shù)，求最經(jīng)濟(jì)的訂購(gòu)方案。設(shè)mij表示第i個(gè)供貨商在第j周的訂貨量，因?yàn)樵撈髽I(yè)對(duì)供應(yīng)商實(shí)際提供的原材料總是全部收購(gòu)，如果第i個(gè)供貨商在第j周被選上，則mij=dij，如果第i個(gè)供貨商在第j周沒(méi)有被選上，則mij=0。根據(jù)前面的決策變量xij的定義，有mij=xij×dij。因?yàn)槿愒牧线\(yùn)輸和儲(chǔ)存的單位費(fèi)用相同，不妨設(shè)運(yùn)輸單位費(fèi)用為a,儲(chǔ)存的單位費(fèi)用為b，假設(shè)a+b與原材料C的單價(jià)比是1:1，若原材料C的單價(jià)為1，原材料A和原材料B的單價(jià)分別為1.2和1.1，則M家供貨商在N周內(nèi)原材料訂購(gòu)的費(fèi)用為：

其中TA，TB，TC是原材料A，B，C三類供貨商集合分別對(duì)應(yīng)的指標(biāo)集。

第三步，以盡量少地采購(gòu)A類和盡量多地采購(gòu)C類為目標(biāo)函數(shù)。通過(guò)增加權(quán)重轉(zhuǎn)化成如下的單目標(biāo)函數(shù)：

2.2 寫(xiě)出約束條件

因?yàn)榭偟墓┴浟繒?huì)在轉(zhuǎn)運(yùn)的過(guò)程中有損耗，根據(jù)近5年8家轉(zhuǎn)運(yùn)商的相關(guān)數(shù)據(jù)，平均的損耗率是1%，所以這里假設(shè)生產(chǎn)企業(yè)的接收量是供貨量的99%。根據(jù)題目的說(shuō)明，為了保障企業(yè)正常生產(chǎn)的需要，該企業(yè)要盡可能保持不少于兩周生產(chǎn)需求的原材料庫(kù)存量，企業(yè)每周產(chǎn)量是28200立方米，兩周是56400立方米。按1%的損耗率得到約束條件為：

2.3 整體模型

記A=(aij)M×N，設(shè)B=[28484.8×2,28484.8×3,…,28484.8×25]T，B中第j個(gè)分量為Bj。再記A=(A1,A2,…,AN),X=(X1,X2,…,XN)，則約束條件表達(dá)式(2)整理為：

綜上，該模型為非線性多目標(biāo)0-1規(guī)劃模型：minf=min(z1,z2,z3)，其中：

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

3.1 V型二進(jìn)制粒子群算法（V-BPSO）

3.1.1 二進(jìn)制粒子群算法（BPSO）

二進(jìn)制粒子群算法（BPSO）最初是在1997年由J.Kennedy和R.C.Eberhart設(shè)計(jì)，是在離散粒子群算法基礎(chǔ)上，約定位置向量、速度向量均由0、1值構(gòu)成。算法流程如下：

初始化粒子位置：按一定策略，生成二進(jìn)制編碼；

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

速度更新公式：

位置的更新公式：

概率映射方式，采用sigmoid函數(shù)將速度映射到[0,1]區(qū)間作為概率，這個(gè)概率就是粒子下一步取值為1的概率；

位置變化的絕對(duì)概率：

當(dāng)前位置為0變化為1，當(dāng)前為1變化為0，這二者被稱為絕對(duì)變化；概率表示為：

公式(6)中，i=1,2,…,M，M為該群體中的粒子總數(shù)，其中ω為非負(fù)的慣性因子，學(xué)習(xí)因子c1和c2是[0,2]之間的常數(shù)；rand1()和rand2()是介于[0,1]之間的隨機(jī)數(shù)；Pbestid是粒子i在第d維的個(gè)體極值點(diǎn)的位置；Gbestd是整個(gè)種群在第d維的全局極值點(diǎn)的位置。

3.1.2 V型二進(jìn)制粒子群算法（VBPSO）

BPSO有很強(qiáng)全局搜索能力，但不能收斂于全局最優(yōu)值，且隨著算法迭代搜索隨機(jī)性越來(lái)越強(qiáng)，缺乏后期的局部搜索能力，因此有學(xué)者對(duì)BPSO算法進(jìn)一步改進(jìn)，Seyedali Mirjalili，Andrew Lewis提出二元粒子的S型和V型傳遞函數(shù)[5]，通過(guò)引入一個(gè)新的傳遞函數(shù)族來(lái)改進(jìn)BPSO算法，在本文中，采用的位置更新函數(shù)公式為：

這里的(t)和(t)表示粒子i在第k維迭代t處的位置和速度，((t))-1是(t)的補(bǔ)碼，這種方法的優(yōu)點(diǎn)是這個(gè)傳遞函數(shù)不會(huì)強(qiáng)制粒子取0或1的值，它們鼓勵(lì)粒子在速度值較低時(shí)停留在當(dāng)前位置，或者在速度值較高時(shí)切換到它們的補(bǔ)碼，（9）式被稱為V型傳遞函數(shù)，粒子的位置更新圖如圖1所示。

圖1 粒子的位置更新

3.2 非支配排序遺傳算法（NSGA-II）

非支配排序遺傳算法[6]是解決多目標(biāo)優(yōu)化的一種算法，它在選擇、交叉、變異算子與普通的遺傳算法相同，在這基礎(chǔ)上，非支配排序遺傳算法增加了兩個(gè)部分，一個(gè)是非支配排序，另一個(gè)是定義擁擠距離及擁擠比較算子。

（1）非支配排序

在非支配排序算法中，每個(gè)解都會(huì)被分配兩個(gè)參數(shù)np和Sp，其中np表示支配第p個(gè)解的個(gè)數(shù)，Sp為種群中被第p個(gè)解支配的個(gè)體集合。非支配排序算法的主要步驟為：

Step1：對(duì)于種群中np=0的個(gè)體，保存在當(dāng)前集合F1中，并記其支配等級(jí)為1；

Step2：對(duì)于F1中的個(gè)體i，其所支配的個(gè)體集合為Si，遍歷Si中的每個(gè)個(gè)體p，執(zhí)行np=np-1，如果np=0則將個(gè)體i保存在集合H中，并記其支配等級(jí)為2；

Step3：以H作為當(dāng)前集合，重復(fù)上述步驟，找到所有的非支配層。

（2）擁擠距離

擁擠距離描述群體個(gè)體擁擠程度（密度）和目標(biāo)函數(shù)值有關(guān)。計(jì)算方法如下：

首先設(shè)個(gè)體的擁擠距離用d表示，個(gè)體i的擁擠距離用di表示，設(shè)置di=0，再設(shè)fm為目標(biāo)函數(shù)，m=1,2,..M，將fm升序排列，第一位和最后一位的擁擠距離設(shè)為d1=d L=∞，對(duì)于非邊界個(gè)體i擁擠距離度量值的計(jì)算公式為：

（3）擁擠比較算子

通過(guò)非支配排序以及擁擠度計(jì)算之后，種群中的每個(gè)個(gè)體n都得到兩個(gè)屬性：非支配序nrank和擁擠度nd。定義擁擠度比較算子為≥,個(gè)體優(yōu)劣的比較依據(jù)為：i≥j，即個(gè)體i優(yōu)于個(gè)體j，當(dāng)且僅當(dāng)irank≤jrank且id≥jd。

4 模型的求解及Pareto前沿圖

本次求解實(shí)驗(yàn)基于MATLAB-R2020b進(jìn)行編程，運(yùn)行環(huán)境為：Windows 10(64位操作系統(tǒng))，內(nèi)存16GB處理器為Intel(R)Core(TM)i7-10510U CPU@1.80GHz 2.30 GHz。

首先設(shè)置VBPSO算法的參數(shù)，種群規(guī)模Np=150，最大迭代次數(shù)maxgen=500，慣性系數(shù)的初始值c1=c2=w=2，速度最大值vmax=6，最大慣性權(quán)重ωmax=0.9，最小慣性權(quán)重ωmin=0.4，其次對(duì)于約束條件利用外點(diǎn)法設(shè)置懲罰函數(shù)。實(shí)驗(yàn)數(shù)據(jù)見(jiàn)參考文獻(xiàn)[7]，運(yùn)行算法所得到Pareto前沿圖如圖2、3所示。

圖2 VBPSO Prareto前沿

圖3 NSGA-IIPrareto前沿

通過(guò)觀察Pareto最優(yōu)解可知，f2和f3兩個(gè)優(yōu)化目標(biāo)之間存在著明顯的約束關(guān)系，當(dāng)f2減小時(shí)，f3是逐漸增大的，從而可證明本文所選目標(biāo)函數(shù)符合多目標(biāo)優(yōu)化目標(biāo)函數(shù)要求。這兩種算法相比，VBPSO運(yùn)行的速度比NSGA-II快，但是后者的Pareto前沿的點(diǎn)分布更均勻，計(jì)算更準(zhǔn)確。

5 用R方法選擇最佳解

在求得多目標(biāo)優(yōu)化的Pareto最優(yōu)解后，決策者可以從可用的Pareto最優(yōu)解方案中找到最佳替代解決方案。對(duì)于每一個(gè)Pareto最優(yōu)解，優(yōu)化目標(biāo)之間的折中是不同的，這里用R方法尋找這種最佳的折中解決方案，具體的步驟如下：

步驟1：列一個(gè)決策表，包含與目標(biāo)相對(duì)應(yīng)的Pareto最優(yōu)解方案的性能數(shù)據(jù)。

步驟2：根據(jù)決策者認(rèn)為的1，2，3等的重要性對(duì)目標(biāo)進(jìn)行排名。如果兩個(gè)或多個(gè)目標(biāo)被認(rèn)為同等重要，則為這些目標(biāo)分配平均排名。

步驟3：根據(jù)目標(biāo)相關(guān)的性能數(shù)據(jù)，按1，2，3等對(duì)Pareto最優(yōu)解方案進(jìn)行排名。如果兩個(gè)或者多個(gè)方案具有相應(yīng)于目標(biāo)的相同值，則將平均排名分配給這些方案。

步驟4：將分配給目標(biāo)等級(jí)和Pareto最優(yōu)解方案轉(zhuǎn)換為相應(yīng)的權(quán)重。轉(zhuǎn)換的公式如下：

wj為目標(biāo)/備選的權(quán)重，rk為目標(biāo)/備選的排序，n為目標(biāo)/備選的數(shù)量。

步驟5：通過(guò)將目標(biāo)權(quán)重與Pareto最優(yōu)解方案的相應(yīng)權(quán)重相加，計(jì)算各方案綜合得分。

步驟6：根據(jù)綜合得分，具有最高綜合得分的方案被認(rèn)為是最佳妥協(xié)Pareto最優(yōu)解。

本模型由VBPSO求得Pareto最優(yōu)解，共15個(gè)備選訂購(gòu)方案，對(duì)每個(gè)方案的三個(gè)目標(biāo)函數(shù)值按數(shù)值的大小進(jìn)行排序，因?yàn)閒1數(shù)值相等，故只需要對(duì)f2和f3排序，同時(shí)按重要性排序f2＞f3＞f1，利用公式（12）轉(zhuǎn)化成相應(yīng)的權(quán)重。

通過(guò)目標(biāo)權(quán)重與備選方案的相應(yīng)權(quán)重相乘后再求和，得到備選方案的綜合得分。其中方案4的綜合得分最高，認(rèn)為它是最佳折中的訂購(gòu)方案，綜合得分見(jiàn)表1。

表1 各方案綜合得分

6 結(jié)語(yǔ)

本文基于生產(chǎn)企業(yè)訂購(gòu)原材料的三個(gè)需求目標(biāo)和該企業(yè)對(duì)供應(yīng)商實(shí)際提供的原材料總是全部收購(gòu)的特點(diǎn)，建立原材料訂購(gòu)方案的模型，最后轉(zhuǎn)化成非線性的0-1規(guī)劃問(wèn)題。該模型采用V型的二進(jìn)制粒子群算法（V-BPSO），引入一個(gè)新的函數(shù)更新粒子的位置，改進(jìn)粒子群算法容易陷入局部最小的缺點(diǎn)，同時(shí)也提高收斂速度。用該算法對(duì)原材料訂購(gòu)模型求解，找到Pareto前沿，用R方法對(duì)備用方案評(píng)價(jià)，得到最佳折中解。因此該模型具有一定的理論價(jià)值和實(shí)踐價(jià)值。