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

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

0 引言

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

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

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

1 作業(yè)調(diào)度模型的建立

1.1 調(diào)度模式分析

當(dāng)前，農(nóng)機服務(wù)公司或農(nóng)機合作社擁有的無人機數(shù)量多，是無人機植保社會化服務(wù)的主力軍[17]。但受成員文化水平及技術(shù)的限制，無人機作業(yè)管理智能化水平較低，傳統(tǒng)的人工調(diào)度經(jīng)驗難于滿足復(fù)雜的植保作業(yè)需求。

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

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

1.2 相關(guān)影響因素數(shù)學(xué)形式表達

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

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

(1)

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

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

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

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

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

1.3 作業(yè)調(diào)度數(shù)學(xué)模型

1.3.1目標函數(shù)

(1) 作業(yè)總收益最大化

(2)

式中F——所有農(nóng)機作業(yè)的總收益，元

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

(2)作業(yè)總時長最小化

minT=maxteij-mintsij-Tr

(3)

式中T——作業(yè)總時長，h

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

1.3.2約束條件

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

(4)

(5)

(6)

tsij

(7)

(8)

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

2 作業(yè)調(diào)度算法設(shè)計

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

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

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

2.1.1影響排序的相關(guān)因素

(1)病蟲害侵染狀況

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

(2)時間窗

農(nóng)戶可根據(jù)以往病蟲害爆發(fā)規(guī)律、當(dāng)年氣候情況、病蟲害擴散趨勢和速度等，設(shè)置訂單作業(yè)時間窗。在訂單優(yōu)先級排序時，訂單要求的計劃開始作業(yè)時間越早，時間窗長度越短，訂單的優(yōu)先級越高。

(3)作業(yè)面積

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

2.1.2算法設(shè)計

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

p=w1Aj+w2Tsj+w3bj

(9)

其中

bj=Tej-Tsj

(10)

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

w1——變量Aj在優(yōu)先級函數(shù)中的權(quán)重

w2——變量Tsj在優(yōu)先級函數(shù)中的權(quán)重

w3——變量bj在優(yōu)先級函數(shù)中的權(quán)重

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

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

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

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

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

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

2.2 基于NSGA-Ⅱ的作業(yè)路徑規(guī)劃算法

2.2.1染色體編碼

染色體編碼是遺傳算法成功實施優(yōu)化的關(guān)鍵。結(jié)合上文排序算法的設(shè)計，此處采用雙層編碼的編碼方式，如圖1所示。

圖1 編碼方式Fig.1 Coding mode

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

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

傳統(tǒng)的遺傳算法中的初始種群多通過隨機生成的方法構(gòu)建，對于飛防隊調(diào)度問題可能會產(chǎn)生大量的劣質(zhì)解或非法解。為避免無效解的生成以及提高算法的運算效率，本文采取以下方法產(chǎn)生初始種群。

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

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

(11)

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

(1)導(dǎo)入完成排序的農(nóng)田訂單集O，以及飛防隊集M，設(shè)置i=j=1。

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

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

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

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

2.2.3模型約束的處理

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

(12)

2.2.4遺傳算子

(1)選擇算子

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

(2)交叉算子

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

若對本文所設(shè)計的編碼結(jié)構(gòu)進行常規(guī)交叉操作，必然會出現(xiàn)大量非法解，增加算法的復(fù)雜度，因此本文采用變異的部分映射交叉算子(Partially mapped crossover，PMX)對所選擇的父代個體進行交叉操作，以提高可行解的比例?；静僮鳛?在一層農(nóng)田編碼中隨機產(chǎn)生兩個交叉點，即交叉點 1 和 2，將兩父代個體A1、A2的對應(yīng)交叉點下層基因?qū)?yīng)交換，生成新個體B1、B2。以含4個農(nóng)田以及4個飛防隊編碼的染色體為例，F(xiàn)2、F4農(nóng)田下的飛防隊發(fā)生交叉，交叉過程如圖2所示。

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

(3)變異算子

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

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

2.2.5整體算法步驟

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

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

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

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

(4)若滿足i>gen，則輸出當(dāng)前最優(yōu)Paroto解集Xi,繪制相關(guān)圖表，結(jié)束程序；否則跳轉(zhuǎn)至步驟(5)。

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

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

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

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

(9)記錄當(dāng)前種群Pi全局最優(yōu)Paroto解集為Xi，i=i+1，跳轉(zhuǎn)至步驟(4)。

3 實例分析

3.1 算例

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

3.2 算法運行結(jié)果分析

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

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

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

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

表1 農(nóng)田作業(yè)訂單信息Tab.1 Information for each job order

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

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

表3 其他相關(guān)作業(yè)信息Tab.3 Other relevant working information

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

表4 Pareto最優(yōu)解集目標值Tab.4 Pareto optimal solution set

圖4 調(diào)度甘特圖Fig.4 Scheduling Gantt chart

圖5 無人機飛防隊轉(zhuǎn)移路線Fig.5 Roadmap of UAV teams transfer

在實際應(yīng)用中，決策者可針對生產(chǎn)需要，從Pareto最優(yōu)解集中選擇最佳作業(yè)方案。

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

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

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

圖7 調(diào)度相關(guān)成本隨迭代過程的變化曲線Fig.7 Changing curve of scheduling related cost with iteration process

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

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

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

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

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

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

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

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

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

因此，在實際生產(chǎn)中，飛防隊可鼓勵無急迫作業(yè)需求的農(nóng)戶，設(shè)置較大的訂單時間窗長度，以降低作業(yè)總成本和作業(yè)總時長，實現(xiàn)更高效合理的調(diào)度。

4 結(jié)論

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

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

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