基于多目標優化的飛防隊作業調度模型研究

曹光喬 張慶凱, 陳 聰 張 萌 張進龍 黃玉祥

(1.農業農村部南京農業機械化研究所， 南京 210014； 2.西北農林科技大學機械與電子工程學院， 陜西楊凌 712100)

0 引言

近年來，無人機植保技術憑借其作業效率高、成本低等特點，發展極為迅速，國內植保無人機保有量及作業面積逐年大幅增長[1-5]。但是，農機社會化組織服務能力弱，農機主管部門對植保飛防隊缺乏有效調度手段，限制了無人機植保技術的進一步推廣[6]。鑒于病蟲害防治的強時效性以及有限的無人機資源，建立高效的智能調度模式對于提高病蟲害防治效率、降低作業成本具有重要意義。

目前，國內外對于農業領域的無人機調度研究相對較少，現有研究多集中于農田內的航線規劃[7-9]，相關農機調度的研究多集中于作物收獲環節。在不同作業場景下，各調度模型所考慮的變量及約束條件是其主要差異。張璠等[10-11]建立了適用于機主選擇的農機調配模型，以作業收益最大作為優化目標，應用基于啟發式優先級規則的農機調配算法進行求解，隨后針對農機應急調度場景，提出了兩種基于優先策略的緊急調配算法，對以作業損失和調度成本為優化目標的調度模型進行了解算。吳才聰等[12]提出帶時間窗的多目標農機調度模型，并利用動態規劃的思想對模型進行了求解。EDWARDS等[13]提出了考慮農田作業條件的農機調度模型，并利用改進的禁忌搜索算法進行求解，該模型適用于多作業環節順序執行的場景。THUANKAEWSING等[14]以產量最高為優化目標，提出了甘蔗收獲機調度模型，該模型以各農田產量最高時段的收獲比例作為約束。HE等[15]建立了以作業總時間為優化目標的調度模型，將各收獲機收獲時間差異作為約束，且考慮了不同機型、農田土壤類型對于調度的影響。與農作物收獲等作業環節相比，植保作業周期更短，且農田病蟲害爆發具有隨機性，各個農田侵染狀況也存在差異。因此，現有農機調度模型及算法不能直接應用于飛防隊作業調度問題。現有研究表明，農機資源調配屬于多目標優化問題[12,16]，但目前研究多將農機資源調度問題轉化為單目標優化問題進行求解[10, 12]。

本文以面向訂單的多飛防隊協同作業模式為對象，對植保飛防隊調度模式進行分析，建立多目標飛防隊作業調度模型，并提出訂單優先級排序算法和基于NSGA-Ⅱ的作業路徑規劃算法，旨在提出合理的飛防隊調配方案，以提高農機的利用效率和作業效益。

1 作業調度模型的建立

1.1 調度模式分析

當前，農機服務公司或農機合作社擁有的無人機數量多，是無人機植保社會化服務的主力軍[17]。但受成員文化水平及技術的限制，無人機作業管理智能化水平較低，傳統的人工調度經驗難于滿足復雜的植保作業需求。

為適應訂單式、托管式、統防統治等農田作業模式，植保飛防隊的作業調度需要解決多個作業訂單的作業排序、作業時間安排、飛防隊作業路徑規劃等問題。飛防隊調度問題的關鍵在于無人機資源與訂單信息的協同整合，生成最優的作業方案。在病蟲害防治作業季，農戶通常需預先向農機服務公司或合作社提交訂單信息，通常包括作業時間窗、作業面積信息、作業價格、作業位置、農田病蟲害侵染情況；調度中心通過管理端，對植保無人機位置信息、種類、作業效率等信息進行匯總；農機服務公司或合作社根據匯總的訂單信息，參照歷史調度經驗和策略制定出調配方案，并組織飛防隊按調配方案的轉移路徑和時間實施作業。飛防隊作業調度問題屬于多目標優化問題，在滿足作業質量情況下，植保服務收益是各經營主體優先考慮的目標，同時更短時間完成所有訂單，可以降低病蟲害造成的作物產量和品質損失。因此，本文的調度目標為在滿足各項約束的情況下總收益最大、作業總時長最小。

結合農業生產需求以及復雜的農田環境，本文研究基于以下假設：①在作業過程中，作業質量不隨調度方案發生變化。②一個飛防隊擁有一臺植保無人機，若干工作人員，由車輛運載轉移。③無人機作業效率恒定且作業過程中無故障。④飛防隊每日可工作時長固定，開展防治作業和田間轉移，其余時間為非作業時間。⑤一個飛防隊可響應多個訂單，一個訂單亦可由多個飛防隊完成，多個飛防隊協同作業時，到達目標農田的時間可不同，當全部作業完成后可同時離開。⑥單個飛防隊的轉移路徑以初始合作社位置為起點，以最后完成的訂單位置為終點。⑦飛防隊在訂單時間窗前到達目標農田，需等待至作業訂單時間窗下限方能開始作業，作業結束時間不得超過訂單時間窗上限。

1.2 相關影響因素數學形式表達

(1)定義集合M={m1,m2,…,mm}表示m組無人機飛防隊，每個飛防隊可以表示為

mi={{xi,yi},w,v,C}

(1)

其中{xi,yi}為飛防隊i的經緯度位置信息；w為飛防隊生產率；v為飛防隊的轉移速度。C={cs,cw,ct,cd}表示飛防隊作業中所產生的各項收入和成本。cs為植保無人機單位面積作業收費；cw為植保無人機單位面積的使用成本，包含機具折舊、飛防隊人工費用和動力費用等；ct為飛防隊單位距離的轉移成本，主要為車輛燃油消耗和駕駛員人工費用；cd為飛防隊單位等待時間的成本。

(2)定義集合O={o1,o2,…,on}表示n個農田訂單作業信息集，單個訂單信息為oj={{xj,yj},Aj,{Tsj,Tej},lj}。其中{xj,yj}為訂單j的經緯度位置信息；Aj為訂單j的農田作業面積；{Tsj,Tej}為訂單j的時間窗，Tsj為計劃開始作業時間，Tej為最晚結束作業時間；lj為農田侵染狀況，其影響訂單優先級順序。

(3)定義集合Pr={Vp,dgh}表示路網中各路徑節點及路徑信息，其中Vp=Vf∪Vm。Vf={Vf1,Vf2,…,Vfm}為農田節點集合；Vm={Vm1,Vm2,…,Vmn} 為各飛防隊初始位置節點集合；dgh為節點g、h之間的距離。

(4)定義集合Ma={pi(g,h)，xij}為相關標志位符號。其中pi(g,h)為路徑轉移標志位,g∈Vf∪Vm,h∈Vf∪Vm, 若飛防隊i經過節點g到達節點h則pi(g,h)為1，否則為0;xij為作業標志位，若飛防隊i在農田j作業則xij為1，否則為0。

(5)定義Tz={ti(g,h),taij,tsij,teij,Tr}為相關時間的集合。其中ti(g,h)為飛防隊i在點g、h之間轉移所消耗的時間；taij為飛防隊i到達農田j的時間；tsij為飛防隊i在農田j的實際開始作業時間；teij為飛防隊i在農田j的實際完成作業時間；Tr為整個作業進程中的非作業時間總和。

1.3 作業調度數學模型

1.3.1目標函數

(1) 作業總收益最大化

(2)

式中F——所有農機作業的總收益，元

農機作業的總收益等于作業的總收入減去總成本。作業總成本包括無人機使用成本、飛防隊等待時間成本以及飛防隊轉移成本。其中無人機使用總成本等于單位面積的使用成本與總作業面積的乘積，等待時間總成本等于單位等待時間成本與所有飛防隊各任務總等待時長的乘積，轉移成本等于所有飛防隊轉移總路程與單位路程轉移成本的乘積。其中等待時間成本及轉移成本為調度相關成本，通過合理的調度方案可降低成本。

(2)作業總時長最小化

minT=maxteij-mintsij-Tr

(3)

式中T——作業總時長，h

作業總時長為最早一個訂單的實際開始作業時間與最后一個訂單的實際完成作業時間之間的時長，減去整個作業進程中總的非作業時間。

1.3.2約束條件

通過對無人機飛防隊調度過程的分析，確定主要約束條件為

(4)

(5)

(6)

tsij

(7)

(8)

式(4)表示所有訂單均有飛防隊進行服務；式(5)～(7)為訂單作業時間的相關約束，式(5)表示飛防隊i經作業點g到達點h的時間等于離開點g的時間加上路程轉移時間；式(6)表示飛防隊i在訂單g的實際完成作業時間等于飛防隊i在點g的實際開始作業時間加上在點g作業的時間；式(7)表示訂單的硬時間窗約束，訂單的實際完成時間不得晚于訂單時間窗要求；式(8)表示進入農田的飛防隊和離開農田的飛防隊數目相等。

2 作業調度算法設計

飛防隊作業調度算法應能為各農田訂單分配合適的飛防隊，同時為飛防隊規劃合理的轉移路線，該調度方案需同時滿足總收益最大、作業總時長最小兩個優化目標，因此飛防隊作業調度問題屬于多目標優化問題。大多數情況下多目標優化問題不存在同時滿足所有目標最優的解，各優化目標之間會相互沖突，只能協調各優化目標，最優解并不唯一,而是Paroto解集，需由決策者進行均衡[18-20]。針對此類復雜的優化問題，傳統的方法如線性加權法、約束法等往往將多目標轉化為單目標進行處理，但目標權重難以確定[21-22]。目前用于求解多目標優化問題Paroto解集的算法有：遺傳算法[23-24]、禁忌搜索算法[25-26]、粒子群算法[27]、蟻群算法[28-29]等。其中Deb提出的帶精英策略的非支配性排序的遺傳算法(NSGA-Ⅱ)因具有良好的分布性和較快的收斂速度，被廣泛應用于各類優化問題分析。本文考慮病蟲害防治需求及算法求解效率，設計了考慮病蟲害程度的作業排序算法和基于NSGA-Ⅱ的作業路徑規劃算法，分兩步對飛防隊作業調度模型進行求解。

2.1 考慮病蟲害侵染狀況的訂單優先級排序算法

病蟲害防治具有強時效性，本文所設計的目標函數以及約束也與各訂單的作業順序緊密相關。在調度時，按照一定的優先級規則對訂單進行排序，然后按次序進行調度作業，可獲得較好的優化目標函數值，且更適應病蟲害防治的實際需要。

2.1.1影響排序的相關因素

(1)病蟲害侵染狀況

現有研究較少涉及病蟲害防治的適時性損失，較難將病蟲害的爆發風險或經濟損失引入調度模型，但經驗表明對病蟲害嚴重的訂單優先作業，可有效降低病蟲害擴散風險及經濟損失。因此，本文將病蟲害侵染狀況作為訂單作業排序的關鍵因素。農戶在提交訂單時，可按照常規觀測方法，將作業訂單的病蟲害等級，設置為重度、中度和輕度3個等級。

(2)時間窗

農戶可根據以往病蟲害爆發規律、當年氣候情況、病蟲害擴散趨勢和速度等，設置訂單作業時間窗。在訂單優先級排序時，訂單要求的計劃開始作業時間越早，時間窗長度越短，訂單的優先級越高。

(3)作業面積

農業病蟲害具有擴散性，對連片面積較大的訂單優先作業，可以獲得良好的防治效果，亦更能發揮無人機高效作業的優勢。因此，農戶訂單作業面積越大，訂單優先級越高。

2.1.2算法設計

本文所設計的訂單作業排序算法分為組間排序與組內排序兩個步驟。首先按照農田侵染狀況對訂單進行分組，然后按照訂單時間窗和作業面積的優先級函數進行組內排序，該優先級函數為

p=w1Aj+w2Tsj+w3bj

(9)

其中

bj=Tej-Tsj

(10)

式中bj——訂單j的時間窗長度，d

w1——變量Aj在優先級函數中的權重

w2——變量Tsj在優先級函數中的權重

w3——變量bj在優先級函數中的權重

在計算時考慮到量綱的統一，式(9)中Aj、Tsj、bj均進行歸一化處理，并有w1+w2+w3=1，各變量權重可根據實際調度需求進行調整。

綜上所述，作業排序算法步驟如下：

(1)取出農田作業訂單集O，并建立3個訂單類別子集O1、O2、O3。

(2)將所有設置為重度病蟲害的農田訂單添加至O1，設置為中度病蟲害的農田訂單添加至O2，其余訂單添加至O3。

(3)在O1、O2、O3內分別計算每個訂單ok的優先值pk，并按照pk大小將訂單排序。

(4)將O1、O2、O3依次連接，對所有訂單重新編號，完成農田作業訂單集O優先級排序。

2.2 基于NSGA-Ⅱ的作業路徑規劃算法

2.2.1染色體編碼

染色體編碼是遺傳算法成功實施優化的關鍵。結合上文排序算法的設計，此處采用雙層編碼的編碼方式，如圖1所示。

圖1 編碼方式Fig.1 Coding mode

第1層為農田編碼，第2層為飛防隊編碼，圖1中第1層編碼為農田F1～F4，第2層中各元素代表在上層農田中作業的飛防隊編號，如在F1訂單中作業的飛防隊為m1、m2，同一訂單中的飛防隊編號不分先后順序。

2.2.2基于貪婪思想的種群初始化

傳統的遺傳算法中的初始種群多通過隨機生成的方法構建，對于飛防隊調度問題可能會產生大量的劣質解或非法解。為避免無效解的生成以及提高算法的運算效率，本文采取以下方法產生初始種群。

定義：設Tc(i,j)為飛防隊i到達農田j時的適時度，Tc(i,j)越小，適時度越高，其中

Tc(i,j)=|Tsj-taij|

(11)

為農田選擇適時度高的飛防隊，可獲得較好的兩目標函數值，種群初始化流程如下：

(1)導入完成排序的農田訂單集O，以及飛防隊集M，設置i=j=1。

(2)判斷O是否為空，若為空則跳轉至步驟(6)，否則順序執行。

(3)依序取出訂單oj，按照公式計算所有飛防隊相對于訂單oj的Tc(i,j)，并按照Tc(i,j)由小到大的次序將飛防隊m排序。

(4)取出Tc(i,j)值最小的飛防隊mi分配給oj。

(6)整理工作表，結束。

2.2.3模型約束的處理

考慮到本文使用的時間窗為單邊硬時間窗，若直接采用所建模型的目標函數構建適應度函數，在進化過程中會產生較多無效解。為保證種群數量，本文采用罰函數法對時間窗約束進行處理，即通過給超出訂單最晚作業時間的解的目標函數一定的懲罰值，以降低不可行解進行遺傳操作的概率，可達到降低不可行解在種群中的比例[30]。在算法中需對優化目標添加懲罰項，式(3)則變為

(12)

2.2.4遺傳算子

(1)選擇算子

普通遺傳算法可利用目標函數所構造的適應度函數進行父代種群的選取，NSGA-Ⅱ從第2代起則需將子代種群與父代種群合并，根據各染色體之間的非支配性關系，同一層級個體之間的擁擠度，選取父代種群進行雜交和變異操作。本文使用標準的快速非支配性排序算法和擁擠度計算算法。

(2)交叉算子

交叉算子影響遺傳算法在解空間中的搜索能力，對于遺傳算法達到全局最優起著關鍵作用[31]。

若對本文所設計的編碼結構進行常規交叉操作，必然會出現大量非法解，增加算法的復雜度，因此本文采用變異的部分映射交叉算子(Partially mapped crossover，PMX)對所選擇的父代個體進行交叉操作，以提高可行解的比例。基本操作為:在一層農田編碼中隨機產生兩個交叉點，即交叉點 1 和 2，將兩父代個體A1、A2的對應交叉點下層基因對應交換，生成新個體B1、B2。以含4個農田以及4個飛防隊編碼的染色體為例，F2、F4農田下的飛防隊發生交叉，交叉過程如圖2所示。

圖2 染色體交叉過程Fig.2 Chromosome crossing process

(3)變異算子

變異操作是保證種群多樣性的重要手段，本文變異算子的基本操作為:在選中個體農田編碼中隨機產生1個變異點，對該農田編碼下層農機編碼片段進行隨機變異，生成新個體，規定至少有一個農機編碼片段發生變異，同樣以含4個農田以及4個飛防隊編碼的染色體為例，F3農田中的m3號農機位發生變異，變異過程如圖3所示。

圖3 基因變異過程Fig.3 Gene mutation process

2.2.5整體算法步驟

基于NSGA-Ⅱ算法的飛防隊調度算法運行流程如下：

(1)讀取完成排序的農田作業訂單集O及飛防隊集M，初始化路徑節點實際行車距離矩陣D。

(2) 設置最大迭代次數gen，種群規模Size，交叉概率Pc，變異概率Pm，迭代次數i=1。

(3)按照初始解生成算法構建初始解，隨機變異產生初始種群P0=(x1，x2，…,xn)。

(4)若滿足i>gen，則輸出當前最優Paroto解集Xi,繪制相關圖表，結束程序；否則跳轉至步驟(5)。

(5)按照交叉概率Pc選取進行交叉的父代個體，使用PMX算子對父代個體進行交叉操作。

(6)按照變異概率Pm選取進行變異的個體，對其中的農機編碼片段進行變異操作,得到子代種群Si。

(7)將父代種群Pi-1與子代種群Si合并，構成新的規模為2Size的種群ConX，計算ConX中個體對應的作業總收益和作業總時間，并進行非支配性排序，并計算同支配序列擁擠度。

(8)將種群規模恢復為Size。若ConX中非支配序列為1的個體數大于種群規模Size，則選擇擁擠度較小的個體進入新的父代種群Pi；若數量小于Size，將ConX中非支配序列為1的個體復制到Pi；對非支配序列為2及以上的個體隨機進行兩兩比較，根據非支配等級和擁擠度選取進入種群Pi的個體，直至種群規模恢復為Size。

(9)記錄當前種群Pi全局最優Paroto解集為Xi，i=i+1，跳轉至步驟(4)。

3 實例分析

3.1 算例

以陜西省武功縣小麥“一噴三防”作業為例進行實驗驗證。該項作業需要調度武功縣及周邊乾縣和興平縣等3個合作社的15支飛防隊，為21個基層村的小麥植保提供統防統治作業，作業時限為7 d。基層村歷史訂單信息通過陜西省某植保無人機調度中心信息管理平臺獲取，飛防隊實際作業能力及各項成本信息通過調研合作社獲取，各路徑節點信息及實際行車距離由天地圖API獲取。具體作業訂單信息、無人機植保飛防隊信息、其他要素信息等如表1～3所示。

3.2 算法運行結果分析

在IntelCorei5 CPU 3.0 GHz、內存為8.0 GB、操作系統為Windows 10的個人計算機上采用Matlab R2018a軟件編程實現本文所設計的作業優先級排序算法及飛防隊作業路徑規劃算法。實驗中設置算法的相關參數為：種群規模200，最大進化代數350，交叉率0.8，變異率0.1。

通過代入實例數據，算法可求出一個Pareto最優解集，能獲得各優化目標下的非支配解。表4為上文算例運行一次后的實驗結果，Pareto最優解集中包括3個調度方案，各個方案中的優化目標值如表4所示，所有訂單均能完成作業，同時滿足單邊時間窗的約束。

圖4為獲得最大作業總收益時的各飛防隊調度甘特圖，圖內條塊表示每個飛防隊的作業計劃，包括作業訂單編號、作業開始時間、作業結束時間。

通過本文所設計的算法可得到每個合作社每個飛防隊的轉移路線。以2號農機合作社為例，圖5為4個飛防隊的作業轉移路線圖，其中農機合作社及農田位置由經緯度坐標點給出。 由圖可直觀看出，每個無人機所分配到的作業任務之間具有空間上的鄰近性，表明了實驗結果的合理性。

表1 農田作業訂單信息Tab.1 Information for each job order

注：感染情況隨機生成，0表示為輕度病蟲害，1為中度病蟲害。

表2 合作社飛防隊信息Tab.2 UAV plant protection team informations of cooperatives

表3 其他相關作業信息Tab.3 Other relevant working information

注：此處作業成本僅包括機器折舊、人工成本等，不含藥劑費用。

表4 Pareto最優解集目標值Tab.4 Pareto optimal solution set

圖4 調度甘特圖Fig.4 Scheduling Gantt chart

圖5 無人機飛防隊轉移路線Fig.5 Roadmap of UAV teams transfer

在實際應用中，決策者可針對生產需要，從Pareto最優解集中選擇最佳作業方案。

圖6為Pareto解集在調度相關成本-作業總時長空間的分布情況，非支配解集中的各解相對分散，可知算法可在解空間中實現有效搜索。

圖6 非支配解的分布情況Fig.6 Distribution of non-dominated solutions

運行過程中的各目標函數最優值的變化情況如圖7、8所示，調度相關成本、作業總時長隨迭代次數的增加而降低，且在迭代次數增加至200次左右即可收斂至穩定值。計算結果表明該算法可實現穩定收斂并具有較好的搜索性能。

圖7 調度相關成本隨迭代過程的變化曲線Fig.7 Changing curve of scheduling related cost with iteration process

圖8 作業總時長隨迭代次數的變化曲線Fig.8 Changing curve of total time of operations varied with iterative process

圖9 運算時間隨訂單數的變化曲線Fig.9 Changing curves of operating time of each case

調研數據表明，在單次防治作業任務中，合作社收到村級訂單數量一般不超過20個，單個訂單作業量33.33～200.00 hm2。在實例數據中分別挑選5、10、15、20個訂單數據組成4個算例，運行本文設計的飛防隊調度算法，記錄兩目標函數值優化至穩定值時所花費的時間，如圖9所示，運算時間隨訂單數量的增加而增大，且近似線性增長。當算例訂單數量達到20個時，訂單作業總面積1 976.30 hm2，運算時間為613 s，基本滿足實際調度需求。在進一步的研究中，可結合其他啟發式搜索策略以提高運算效率，同時通過對訂單任務進行合理合并，也可降低問題規模，減少運算時間。

同時，在縣域范圍內，病蟲害防治單次作業合理周期約為7 d左右，訂單作業窗口可由農戶根據實際需要以及提交訂單的先后順序進行排定。為探究時間窗長度的變化對調度結果的影響，本文設計了以下實驗進行驗證。

針對上述實例，保持其他數據不變，分別設置所有訂單時間窗長度為3～5 d，在每個時間窗長度水平下使用本文算法分別進行5次仿真，可得時間窗長度與調度相關成本和作業總時長的關系，仿真結果如表5所示。

表5 時間窗口長度對實驗結果的影響Tab.5 Effect of time window duration on experimental results

注：各項結果為5次實驗的平均值。

由表5可知，在總作業周期不變的情況下，當訂單時間窗口長度大于3 d時，等待時間成本消失，調度相關成本僅為轉移成本，且隨時間窗長度的增加，調度相關成本與作業總時長均呈下降趨勢。

當總作業周期固定時，訂單時間窗口越長，各訂單的時間窗口重疊度增加，時間窗對飛防隊作業調度的約束會降低，不需要為滿足部分訂單的時間窗要求而增加飛防隊的轉移距離；同時訂單時間窗重疊度的增加，降低了因訂單時間窗過于集中或分散造成飛防隊不能滿足訂單作業需求或飛防隊等待作業時間過長的風險，提高了無人機利用效率。

因此，在實際生產中，飛防隊可鼓勵無急迫作業需求的農戶，設置較大的訂單時間窗長度，以降低作業總成本和作業總時長，實現更高效合理的調度。

4 結論

(1)以面向訂單的多飛防隊協同調度為研究對象，對飛防隊作業調度的各項成本進行分析，構建了以調度總收益最大和調度總時長最小為優化目標的飛防隊作業調度模型。該模型考慮了農田面積、無人機和農田位置信息、作業時間窗等因素，滿足多目標調度需求以及單邊硬時間窗的約束，提高了飛防隊作業調度模型的準確性。

(2)通過分析對比多目標調度算法特點，結合飛防隊作業調度需求，提出了考慮病蟲害侵染狀況的作業優先級排序算法和基于NSGA-Ⅱ的作業路徑規劃算法，并對飛防隊作業調度模型進行了求解。

(3)通過Matlab軟件進行了算例及算法運行相關實驗分析，本文提出的飛防隊作業調度模型及算法，求解出滿足時間窗約束的Pareto解集，為決策者提供了多個備選的調度方案。