一種求解資源約束項目調度問題的改進引力搜索算法

劉永利，張曉陽

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

0 引 言

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

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

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

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

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

1 問題描述

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

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

約束條件為

決策變量ajt定義為

RCPSP的目標函數為最小完工時間，即式（1）中的Cmax；式（2）保證每個活動只能被執行一次；式（3）保證活動i在其所有前置活動完成前不能開始執行；式（4）確?；顒颖仨氃谄渌璧目稍偕Y源可用時執行。

2 求解RCPSP的引力搜索算法

2.1 算法簡介

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

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

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

向心力公式為

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

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

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

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

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

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

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

2.2 相關機制

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

2.2.1 調度生成機制

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

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

（1）初始化，n＝1。

（2）計算可調度活動決策集合Dn。

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

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

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

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

（7）n＝n+1。

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

2.2.2 鄰居生成機制

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

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

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

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

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

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

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

2.2.3 混沌機制

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

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

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

圖3 控制參數取值Fig.3 Value of the control parameter

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

2.3 求解步驟

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

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

（2）隨機初始化，并對初始化后的解執行串行調度操作使其可行。

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

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

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

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

3 結果與分析

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

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

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

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

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

表2 實驗結果Tab.2 Experimental results

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

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

3.1 資源利用率分析

資源利用率可評估算法在最小完工時間內的資源利用情況，避免資源浪費。取每個算法運行10次實例結果中的1次，以測試集J3001_01和J12001_01為例，分析PSO，CS，GSA，NCM-IGSA和IGSA在最小完工時間（最優值）內的資源利用率，實驗結果分別如圖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種對比算法得到的最優值分別為45，45，45，43和43。由圖5可以看出，在最優值的時間范圍內，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種算法得到的最優值分別為124，121，120，124和116，同理，由圖6可以看出，測試集為J12001_01時，5種對比算法資源利用率超過0.8的分別有16份，19份，22份，13份和出，PSO在第90次迭代時得到其最優值45，在第228次迭代時相對誤差MPE為0；CS在第91次迭代得到其最優值45，在第266次迭代時相對誤差MPE為0；GSA在第105次迭代得到其最優值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 收斂性分析

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

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

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

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

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

4 結 語

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