劉永利，張曉陽

（河南理工大學 計算機科學與技術學院，河南 焦作 454000）

0 引 言

資源約束項目調(diào)度問題（resource constrained project scheduling problem，RCPSP）是指在可再生與不可再生資源的約束條件下，根據(jù)活動的優(yōu)先關系計算項目的最小完工時間［1］，其在項目管理、施工和生產(chǎn)計劃等諸多領域［2］備受關注。

求解RCPSP主要采用精確算法和啟發(fā)式算法。精確算法包含分支定界［3］、分支裁剪［4］和基于事件的方法［5］。受限于計算復雜度和NP-hard等因素，精確算法僅適用于解決小規(guī)模問題，截至目前，最好的精確算法也只能在實例不受高度資源約束的情況下解決最多涉及60個活動的項目［6］，但現(xiàn)實中項目數(shù)往往超過這個規(guī)模，且要求解決方案必須快速確定，所以針對RCPSP的研究重心集中在啟發(fā)式算法領域，以期在準確性、計算速度、簡化操作和靈活性之間尋找最佳的平衡［6］。R.Zamani［7］提出了一種基于磁鐵 分頻算子的遺傳算法，通過局部搜索提高算法的性能；還提出了基于遺傳算法和隱式枚舉搜索技術的算法［8］，可有效提高算法的計算速度，并易于得到更優(yōu)的解集；S.Elsayed等［9］提出了一種基于雙算子的綜合進化算法，其中的元啟發(fā)式算法自適應選擇，可避免算法陷入局部最小值；JIA Q等［10］提出了一種改進的粒子群優(yōu)化算法，該算法使用排序優(yōu)先技術表示粒子，并基于雙對齊的局部搜索技術改進算法的求解，在解決項目調(diào)度問題庫（a project scheduling problem library，PSPLIB［11］）

中的問題時，該算法的性能優(yōu)于其他PSO變體；A.Thammano等［12］提出了一種基于遺傳算法的綜合方法，該綜合方法將禁忌搜索［13-14］和模擬退火［15］應用于小部分種群中，以減少計算時間；E.Rashedi等［16］提出了引力搜索算法（gravitational search algorithm，GSA），該算法以自然界中常見的萬有引力現(xiàn)象為靈感，將其轉(zhuǎn)化為隨機搜索最優(yōu)解的過程，利用個體間的引力實現(xiàn)群體間的信息交互，從而找到最優(yōu)解。

引力搜索算法因其流程簡單、參數(shù)較少、易于實現(xiàn)等特點，在已有的啟發(fā)式算法中脫穎而出，目前已廣泛應用于不同類型問題的優(yōu)化求解。然而，引力搜索算法也存在不足，如容易陷入局部最優(yōu)和精度不高等，且算法的收斂性易受到迭代次數(shù)的影響，導致計算成本較高。

鑒于此，為了更有效地求解RCPSP，得到項目的最小完工時間，本文將向心力概念引入GSA中，以有效平衡粒子的探索能力與開發(fā)能力，且避免算法陷入局部最優(yōu)；同時，將Logistic混沌映射引入到RCPSP求解，通過Logistic映射得到的混沌序列具有隨機性和不可預測性，這兩種性質(zhì)可以增加算法找到最優(yōu)解的概率，避免陷入局部最優(yōu)解。

1 問題描述

RCPSP中，假定所有活動和優(yōu)先約束均已知，且必須實現(xiàn)所有活動［17］，每個活動的執(zhí)行均須遵循優(yōu)先關系與資源約束，其中優(yōu)先關系是指活動j必須要等到其所有的前置活動pred（j）結(jié)束之后才會被執(zhí)行。

RCPSP的基本描述如下：項目由非循環(huán)節(jié)點活動（activity-on-node，AON）網(wǎng)絡拓撲表示。網(wǎng)絡拓撲結(jié)構(gòu)包含（N+2）個活動，其中第一個和最后一個活動分別表示開始和結(jié)束活動，稱為虛擬活動，不占用時間和資源。所有活動共享R種可再生資源，其中第k種可再生資源的可利用總量為Rk（k＝1，2，…，R）。第j個活動的持續(xù)時間為dj，對資源k的需求量為rjk。此外，所有參數(shù)的值均為非負整數(shù)。單模的RCPSP數(shù)學模型為

約束條件為

決策變量ajt定義為

RCPSP的目標函數(shù)為最小完工時間，即式（1）中的Cmax；式（2）保證每個活動只能被執(zhí)行一次；式（3）保證活動i在其所有前置活動完成前不能開始執(zhí)行；式（4）確保活動必須在其所需的可再生資源可用時執(zhí)行。

2 求解RCPSP的引力搜索算法

2.1 算法簡介

引力搜索算法的靈感來源于萬有引力定律和牛頓第二定律。依照萬有引力定律，搜索空間中的各個粒子之間存在相互作用力并在萬有引力的作用下朝著質(zhì)量較大的粒子方向移動。同時，搜索空間中個體的運動遵循牛頓第二定律。隨著算法不斷迭代，最終整個種群會聚集在質(zhì)量最大的個體周圍，找到的質(zhì)量最大個體即為全局最優(yōu)解。因此，根據(jù)引力搜索算法，可將求解最大值優(yōu)化問題轉(zhuǎn)化為在搜索空間中搜索質(zhì)量最大的個體問題。

相較于其他元啟發(fā)式優(yōu)化算法，GSA算法具有較強的全局搜索能力且操作簡單、易于擴展。但GSA算法的收斂性會受到迭代次數(shù)的增加而降低，并有陷入局部最優(yōu)的傾向，而且當涉及到矩陣運算時GSA算法計算成本較高。針對該問題，提出改進的引力搜索算法（improved gravitational search algorithm，IGSA），該算法將向心力概念引入GSA，可有效避免算法陷入局部最優(yōu)并提高求解精度。

IGSA算法中，向心力大小與物體質(zhì)量、線速度成正比，與物體間距離成反比。質(zhì)量越大，行星朝太陽運動的速度越快，能更快地找到最優(yōu)解。同理，線速度越大和距離越小也能達到這種效果。

向心力公式為

式中：F為物體之間的向心力；M為物體質(zhì)量；V為線速度；R為兩者之間的距離。

根據(jù)式（7），得到線速度V表達式，即

為方便計算，將線速度V表達式更新為

根據(jù)式（9），迭代過程中速度更新公式為

式中：fit為當前搜索代理的適應值；fitg為迭代過程中的最優(yōu)適應值；Dis（x，xg）為當前搜索代理與全局最優(yōu)搜索代理之間的距離。

根據(jù)速度表達式，可得到位置更新公式，即

為了更全面描述搜索過程，該算法將尋找最優(yōu)解的過程分為勘探階段，開發(fā)階段和避免陷入局部最優(yōu)解階段。勘探階段為兩個物體在向心力作用下做近似圓周運動；開發(fā)階段為隨著迭代次數(shù)增加，物體之間的距離逐漸縮??；避免陷入局部最優(yōu)解階段物體之間的距離為定值。開發(fā)階段式（11）可提高算法的收斂速度，并避免陷入局部最優(yōu)。

2.2 相關機制

在求解RCPSP過程中，調(diào)度生成機制與鄰居生成機制必不可少。此外，為進一步提高尋優(yōu)的效率和效果，在求解過程中引入混沌機制，增強算法的多樣性。多樣性是指在生成解決方案后算法有一定的概率可重新執(zhí)行插入操作或交換操作，執(zhí)行完此變異操作后，算法可能會得到更好的解決方案和更優(yōu)的最小完工時間。

2.2.1 調(diào)度生成機制

調(diào)度生成機制通過給定的優(yōu)先級列表為每個活動分配開始時間，從而確定可行的調(diào)度。該機制主要包括串行調(diào)度和并行調(diào)度。本文采用每次執(zhí)行一個活動的方式依次調(diào)度，隸屬于串行調(diào)度。

串行調(diào)度方案通過階段式的擴展部分調(diào)度表生成可行調(diào)度方案，其中部分調(diào)度表是指只有一部分活動被分配完成時間的調(diào)度表［18］。在每個階段，調(diào)度生成方案得到所有可調(diào)度活動的集合稱為決策集，然后使用特定的優(yōu)先級規(guī)則，從決策集中選擇一個或多個活動進行調(diào)度，最終得到調(diào)度結(jié)果。串行調(diào)度的具體步驟如下。

（1）初始化，n＝1。

（2）計算可調(diào)度活動決策集合Dn。

（3）選擇Dn中具有最高優(yōu)先權值的活動i，同時獲取該活動i的持續(xù)時間di和在t時間資源r的剩余量（t＝1，2，…，T；r∈R，t為當前時刻，T為總時間）。

（4）從已調(diào)度集合Sn中獲取活動i的所有前置活動中的最遲完成時間，并將其作為i的最早可能開始時間ESi。

（5）計算活動i在滿足時序和資源約束條件時的最早開始時間Si，并計算活動i的完成時間Fi。

（6）將活動i從Dn中刪除，同時更新Sn+1。

（7）n＝n+1。

（8）若n值不大于活動數(shù)，則跳至步驟（2）。

2.2.2 鄰居生成機制

當前解決方案的鄰居可通過插入或交換操作得到。當執(zhí)行完插入或交換操作后，可能會違反優(yōu)先規(guī)則和資源約束，所以需要重新對新的解決方案進行串行調(diào)度操作使其可行。下面以一個具體實例說明插入和交換操作。表1為當前解決方案，圖1和圖2分別表示對當前解決方案執(zhí)行插入和交換操作。

表1 當前解決方案Tab.1 Current solution

假設實例共9個活動，對表1所示的解決方案分別執(zhí)行插入和交換操作。

（1）插入操作：選擇活動5執(zhí)行插入操作，使其移動到新的位置。在活動5的最后一個前置活動2和最早的后繼活動8之間隨機選擇一個活動，以活動8為例，得到此活動的位置為6，即將活動5插入位置6，并依次移動活動8和活動7的位置，最終結(jié)果如圖1所示。

圖1 插入操作示例Fig.1 Insert operation example

（2）交換操作：選擇活動4執(zhí)行交換操作。在活動4的最后一個前置活動2和最早的后繼活動7之間隨機選擇一個活動，以活動3為例，得到此活動的位置為2，即將活動4交換到活動3的位置2，活動3交換到活動4的位置5，其他活動位置不變。最終結(jié)果如圖2所示。

圖2 交換操作示例Fig.2 Exchange operation example

2.2.3 混沌機制

混沌映射用于生成混沌序列，是一種由簡單的確定性系統(tǒng)產(chǎn)生的隨機性序列，具有非線性、隨機性、遍歷性和長期不可預測性等特征。常見的混沌映射有Logistic映射、Gussian映射、Singer映射等。本文選用簡單易用的一維Logistic混沌映射，其數(shù)學表達式為

式中：n為迭代次數(shù)；μ∈［0，4］，為Logistic控制參數(shù)；x∈（0，1）時，Logistic映射處于完全混沌狀態(tài)。

圖3表示初始值x0一定時，隨著控制參數(shù)μ取值變化，迭代過程中可得到的取值范圍。由圖3可以看出，在μ越接近4的地方，x取值范圍越是平均分布在0到1的區(qū)域。圖4表示初始值x0為0.5，迭代次數(shù)為300，控制參數(shù)μ分別為2.5與3.98時x的取值范圍。由圖4可以看出，μ3.98時迭代生成的值處于隨機分布狀態(tài)，μ2.5時，經(jīng)過一定次數(shù)的迭代后，生成的值將收斂到一個特定的數(shù)值。

圖3 控制參數(shù)取值Fig.3 Value of the control parameter

圖4 不同控制參數(shù)迭代后的值Fig.4 Values of different control parameters after iteration

2.3 求解步驟

基于上述3種機制和算法，改進的引力搜索算法解決RCPSP的步驟如下。

（1）初始化種群數(shù)量n、最大迭代數(shù)MaxIt、Logistic控制參數(shù)μ和初始值x0。

（2）隨機初始化，并對初始化后的解執(zhí)行串行調(diào)度操作使其可行。

（3）生成新的解。按照式（10）和（11）得到步長stepsize。若0＜stepsize≤0.3，對當前解執(zhí)行插入操作；若0.3＜stepsize≤1，對當前解執(zhí)行交換操作。生成的新解可能會違反優(yōu)先規(guī)則或資源約束，所以需要對新解執(zhí)行串行調(diào)度操作使其可行。

（4）評估質(zhì)量。根據(jù)優(yōu)先規(guī)則和資源約束得到最佳適應度的值Cmax。若解的適應度值Fj＞Cmax，更新此解的適應度值；若解的適應度值Fj＝Cmax，計算平均資源利用率，將資源利用率更高的保存在迭代中。

（5）多樣化。分別使用隨機函數(shù)與混沌機制中的式（12）生成0～1的參數(shù)K和Pa，若K＞Pa，則對相應的解進行交換操作得到新的解并對新解重新評估質(zhì)量。

（6）當相對誤差MPE的值為0或達到最大迭代數(shù)MaxIt時，輸出Cmax和對應的調(diào)度方案，否則跳到上述（3）。

3 結(jié)果與分析

為了評估IGSA算法求解RCPSP的效果，選擇PSPLIB數(shù)據(jù)集進行對比實驗。PSPLIB數(shù)據(jù)集包含難易程度不同的調(diào)度問題實例，這些實例根據(jù)每個項目所包含的活動數(shù)量分為幾組，包括J30，J60，J90和J120基準測試問題集等。對于該數(shù)據(jù)集，測試集中包含了每個項目的最優(yōu)解，稱之為參考值。對比算法包括粒子群優(yōu)化算法（particle swarm optimization，PSO）、布谷 鳥 搜 索 算 法（cuckoo search，CS）、GSA和無混沌機制的IGSA（NCM-IGSA）。

為了保證實驗結(jié)果的準確性，對所有的對比算法設置種群數(shù)量n為25，最大迭代數(shù)MaxIt為400。PSO，CS算法，GSA和NCM-IGSA中的參數(shù)Pa為固定值0.25，IGSA的Logistic控制參數(shù)μ為4和初始值x0為0.6，PSO算法的速度和位移公式由式（13）～（14）表示，GSA的引力常量G隨著迭代增加而減小，計算如式（15）所示。通過與測試集中的參考值進行比較，使用相對誤差（MPE）、最優(yōu)值（best cost，BC）衡量算法的性能，其中相對誤差MPE的值等于［（算法求得的最優(yōu)解-參考值）／參考值］。除此之外，為了使實驗結(jié)果更嚴謹，將數(shù)據(jù)集中的實例運行10次，得到平均值（AVG）并進行比較。選擇的測試集為PSPLIB中的J3001_01，J6001_01，J9001_01，J12001_01，其參考值分別為43，77，73和105。實驗結(jié)果如表2所示。

其中，式（13）和（14）分別為在t時刻d維空間中第i個粒子的位置與速度更新公式；w為慣性權重，取w＝1；c1與c2為兩個加速度系數(shù)，用于平衡個人經(jīng)驗和社會經(jīng)驗，一般取c1＝c2＝2；r1和r2為［0，1］內(nèi)均勻分布的隨機數(shù)；pbestd為到目前為止第i個粒子在d維空間中搜索到的最好位置；gbestd為整個粒子群搜索到的全局最優(yōu)位置。

式中：G0為100；α為20；iter為當前迭代次數(shù)；MaxIt為最大迭代數(shù)。

由表2可以看出，PSO，CS，GSA，NCM-IGSA和ICF-GSA 5種算法在測試集J3001_01上皆能達到最優(yōu)值43，相對誤差均為0，但每個算法運行10次后求得的平均值分別為43.6，43.4，43.6，43.4，43.2，由此可見，ICF-GSA算法相較于其他對比算法得到最優(yōu)解的概率更高。所以在解決活動數(shù)為30的小規(guī)模問題時，ICF-GSA的效果優(yōu)于其他4種算法。

表2 實驗結(jié)果Tab.2 Experimental results

在活動數(shù)為60的中規(guī)模調(diào)度問題的測試集J6001_01上，5種算法均能得到最優(yōu)解，相對誤差均為0，但CS，GSA，NCM-IGSA 3種算法的平均值分別為78.8，77.9和78.4，而PSO和ICF-GSA的平均值為0，說明這兩種算法運行10次后的結(jié)果皆為最優(yōu)值77，所以在解決中規(guī)模問題時，PSO和ICF-GSA算法效果優(yōu)于其他3種對比算法。

在解決活動數(shù)為90和120的較大規(guī)模調(diào)度問題測試集J90001_01和J12001_01時，5種算法均未能達到參考值。雖然PSO算法和CS算法得到的最優(yōu)值與IGSA相同，都更接近于參考值，但IGSA的平均值小于PSO算法和CS算法。綜上可知，IGSA在求解RCPSP時較為有效，易于得到較好的調(diào)度結(jié)果。

3.1 資源利用率分析

資源利用率可評估算法在最小完工時間內(nèi)的資源利用情況，避免資源浪費。取每個算法運行10次實例結(jié)果中的1次，以測試集J3001_01和J12001_01為例，分析PSO，CS，GSA，NCM-IGSA和IGSA在最小完工時間（最優(yōu)值）內(nèi)的資源利用率，實驗結(jié)果分別如圖5～6所示，其中最大資源利用率設置為1，紅色實線表示0.8，紫色點劃線表示0.6。

圖5 不同算法在測試集J3001_01上的資源利用率Fig.5 Resource utilization of the different algorithm on the test set J3001_01

當測試集為J3001_01時，5種對比算法得到的最優(yōu)值分別為45，45，45，43和43。由圖5可以看出，在最優(yōu)值的時間范圍內(nèi)，PSO，CS，GSA，NCM-IGSA和IGSA 5種算法資源利用率超過0.6的分別有5份、3份、7份、7份和7份，占總時間的11%，6%，15%，16%和16%，所以NCM-IGSA和IGSA算法能更好地利用資源；當測試集為J120001_01時，5種算法得到的最優(yōu)值分別為124，121，120，124和116，同理，由圖6可以看出，測試集為J12001_01時，5種對比算法資源利用率超過0.8的分別有16份，19份，22份，13份和出，PSO在第90次迭代時得到其最優(yōu)值45，在第228次迭代時相對誤差MPE為0；CS在第91次迭代得到其最優(yōu)值45，在第266次迭代時相對誤差MPE為0；GSA在第105次迭代得到其最優(yōu)值45，在第279次迭代時相對誤差MPE為0；NCM-20份，占總時間的12.9%，15.7%，18.3%，10.5%和17.24%?？梢姡啾绕渌惴?，GSA和IGSA更能充分利用資源，避免造成資源浪費。

圖6 不同算法在測試集J12001_01上的資源利用率Fig.6 Resource utilization of the different algorithm on the test set J12001_01

3.2 收斂性分析

收斂性可衡量算法求解最優(yōu)值的速度。本文同樣以每個算法運行10次實例結(jié)果中的1次且測試集為J3001_01和J6001_01為例，以最優(yōu)值與MPE為度量指標，每種算法迭代400次，分析PSO，CS，GSA，NCM-IGSA和IGSA 5種算法收斂到最優(yōu)解與相對誤差為0的速度。實驗結(jié)果如圖7～8所示。

當測試集為J3001_01時，5種算法得到的最優(yōu)值分別為45，45，45，43和43。由圖7可以看IGSA在第120次迭代得到其最優(yōu)值43，在第294次迭代時相對誤差MPE為0；IGSA經(jīng)過53次迭代得到其最優(yōu)值43，在第209次迭代時相對誤差MPE為0。相較于其他對比算法，IGSA在解決小規(guī)模調(diào)度問題時能更快收斂到最優(yōu)值。

圖7 測試函數(shù)為J3001_01時的收斂曲線Fig.7 Convergence curves with the test function of J3001_01

當測試集為J6001_01時，5種算法得到的最優(yōu)解皆為77，由圖8可以看出，PSO在第36次迭代得到其最優(yōu)值77，在經(jīng)過313次迭代后相對誤差MPE的值為0；CS在第27次迭代得到其最優(yōu)值77，在經(jīng)過223次迭代后相對誤差MPE為0；GSA在第209次迭代得到其最優(yōu)值77，在經(jīng)過367次迭代后相對誤差MPE為0；NCM-IGSA在第108次迭代得到其最優(yōu)值77，在第259次迭代時相對誤差MPE為0；IGSA在迭代剛開始時即能得到最優(yōu)解77，且在第170次迭代時相對誤差MPE為0。相對于其他算法，IGSA在解決中規(guī)模調(diào)度問題時也能更快收斂到最優(yōu)值。這意味著IGSA能夠在搜索空間中快速找到最優(yōu)解，并有效避免陷入局部最優(yōu)，提高了算法的收斂速度。

圖8 測試函數(shù)為J6001_01時的收斂曲線Fig.8 Convergence curves with the test function of J6001_01

4 結(jié) 語

為了提高引力搜索算法的收斂速度并避免陷入局部最優(yōu)，本文在該算法的基礎上引入向心力，并在求解RCPSP的過程中加入混沌機制，提出改進的引力搜索算法IGSA。為了評估算法的優(yōu)化表現(xiàn)，在PSPLIB問題實例J30，J60，J90和J120上進行了實驗。實驗結(jié)果表明，在解決RCPSP問題時，相較于其他算法，IGSA得到的最優(yōu)解更接近于測試集的參考值，在資源利用和收斂性方面也優(yōu)于其他算法。